#   Estudo: Câncer não-melanoma de pele
#    
#   Duas cidades (Mineapolis e Dallas) em 1994 
#    no. de câncer de pele não melanoma
#    por faixa etaria (8)
#     15-24, 25-34, .....,>84
#     controlado pela população (offset)
#=========================================================
#           Entrada e Preliminares dos Dados 
#=========================================================
# 
 contagem<-c(1,16,30,71,102,130,133,40,4,38,119,221,259,310,226,65)
 cidade<-c(rep(0,8),rep(1,8))
 faixa<-c(rep(c(20,30,40,50,60,70,80,90),2))
 pop<-c(172675,123065,96216,92051,72159,54722,32185,8328,181343,
        146207,121374,111353,83004,55932,29007,7538)
#  
#===========================================================
#   Descrição dos Dados 
#=========================================================== 
#
 table(cidade)
 table(cidade,faixa)
#
#=================================================================
#   Modelo de Poisson 
#==================================================================
#
 m1<-glm(contagem~cidade + offset(log(pop)),family=poisson)
 m1
 plot(m1)
#  Incluindo Faixa Etária Contínua
 m2<-glm(contagem~cidade*faixa + offset(log(pop)),family=poisson)
 m2
 plot(m2)
 summary(m2)
#  Incluindo Faixa Etária Categórica
 m3<-glm(contagem~cidade+factor(faixa) + offset(log(pop)),family=poisson)
 plot(m3)
 summary(m3) 
#==========================================================
# Predizer na escala original ---> R group
#      arrumar para o banco de dados
#====================================================================
#
# assuming 'mod' is your fitted model and 'pdata' is the data frame
# containing the two rows of new values you wanted predictions for
 p <- predict(mod2, newdata = pdata, type = "response", se.fit = TRUE)
## apprx 95% CI
 upr <- with(p, fit + (2 * se.fit))
 lwr <- with(p, fit + (2 * se.fit))
## inverse link fun
 invLink <- family(mod)$linkinv
## map these on to the scale of response
 fit <- with(p, invLink(fit))
 upr <- invLink(upr)
 lwr <- invLink(lwr)
============================================================================
           #   Partos Cesarianos
           #    Comparação de Hospitais privados e publicos
           #    Tamanho do hospital (offset)
           #     
           #=========================================================
           #           Entrada e Preliminares dos Dados 
           #=========================================================
           # 
           partos<-c(236,739,970,2371,309,679,26,1272,3246,1904,357,1080,
           1027,28,2507,138,502,1501,2750,192)
           hosp<-c(0,rep(1,5),0,rep(1,6),0,1,0,rep(1,4))
           ces<-c(8,16,15,23,5,13,4,19,33,19,10,16,22,2,22,2,18,21,24,9)
           #  
           pos<-seq(1:20)
           tx<-ces/partos
           plot(tx,main="Taxa de Cesarianas",xlab="Hospital",ylab="Taxa")
           points(pos[hosp==0],tx[hosp==0],pch=16)
           #
           #==========================================================
           #    Explorando Regressão Linear SImples
           #=========================================================== 
           #
           mod1<-lm(ces~partos+hosp) # dois regressores
           plot(mod1)
           summary(mod1) # Absurdo o coef. para hospital é +
           mod2<-lm(sqrt(ces)~partos+hosp) # dois regressores
           plot(mod2)
           summary(mod2) # o coef. para hosp é menor mas continua + 
           mod3<-lm(ces/partos~hosp) # um regressor
           plot(mod3)
           summary(mod3) # agora é negativo
           #
           #=================================================================
           #           Modelo de Poisson 
           #==================================================================
           #
           m1<-glm(ces~hosp+offset(log(partos)),family=poisson)
           summary(m1) # o ajuste náo é adequado mas o coef. é consistente e -
           plot(m1)
           plot(residuals(m1)~fitted(m1))
           abline(h=0) # problema com a homocedasticidade
           #
           plot(log(fitted(m1)),log((ces-fitted(m1))^2),
           xlab=expression(hat(mu)),ylab=expression((y-hat(mu))^2))
           abline(0,1) ## gráfico para verificar se a média = variância
           # Incluindo um parâmetro de dispersão
           (dp<-sum(residuals(m1,type="pearson")^2)/m1$df.res)
           summary(m1,dispersion=dp)
           #  Teste F é mais indicado do que o Qui-quadrado
           drop1(m1,test="F") 
           #  Outra alternativa: variância robusta (estimador sanduiche)
           library(sandwich)
           library(lmtest)
           coeftest(m1, vcov=sandwich)
           ============================================================================
            Modelo Binomial Negativo
           ====================================================================
           require(MASS)
           modn<-glm.nb(ces~offset(log(partos))+hosp)
           summary(modn)
           ===================================================================