CE-225 - Segundo semestre de 2012

No quadro abaixo será anotado o conteúdo dado em cada aula do curso.
É indicado material para leitura correspondente ao conteúdo da aula nas referências bibliográficas básicas do curso:

• T & S: M. Antónia Amaral Turkman & Giovani Loiola Silva
• C & D: Gauss M. Cordeiro & Clarice G. B. Demétrio
• P: Gilberto A. Paula

Considere o conjunto de dados "heart" Sobre número de pacientes com e sem ataque cardíaco para diferentes níveis da enzima creatinina kinase no sangue. O objetivo é examinar se os níveis da enzima podem ser utilizados como auxiliares de diagnóstico.

      Ataque cardíaco
     ------------------
CK		Sim		Nao
-----------------------------------
 20		 2		88
 60		13		26
100		30		 8
140		30		 5
180		21		 0
220		19		 1
260		18		 1
300		13		 1
340		19		 1
380             15		 0
420		 7		 0
460		 8		 0
----------------------------------

Analise os dados. A seguir são dadas algumas (mas não todas!) questões que podem ser consideradas nas análises.

1. Qual modelo (GLM) poderia ser adotado?
2. Que gráficos exploratórios poderiam ser feitos?
3. Como pode-se investigar o efeito da CK ?
4. O modelo NULO é suficiente para explicar os dados?
5. O modelo COMPLETO superior ao NULO para explicar os dados?
6. Como podemos avaliar a qualidade do ajuste do modelo completo?
7. Como podemos detarminar se a covariável é realmente relevante?
8. Como obter intervalos para os parâmetros?
9. Como obter, avaliar e fazer gráficos do modelo predito?
10. Como avaliar diferentes opções para função de ligação?
11. Como avaliar se a linearidade das covariáveis selecionadas? (por exemplo, uma função quadrática ou cúbica no modelo linear seria melhor?)

infarto <- read.table("clipboard")
infarto
names(infarto) <- c("CK","S","N") 
infarto
 
infarto <- transform(infarto, prop = S/(S+N))
 
with(infarto, plot(prop ~ CK))
 
mod0 <- glm(cbind(S,N) ~ 1, family=binomial, data=infarto)
mod0
summary(mod0)
logLik(mod0)
deviance(mod0)
 
mod1 <- glm(cbind(S,N) ~ CK, family=binomial, data=infarto)
mod1
summary(mod1)
logLik(mod1)
deviance(mod1)
anova(mod1)
anova(mod1, test="Chisq")
 
-2*(logLik(mod0) - logLik(mod1))
 
(pred1 <- predict(mod1))
(pl <- with(infarto, coef(mod1)[1] + coef(mod1)[2] * CK))
with(infarto,plot(prop/(1-prop), pred1, asp=1)); abline(0,1)
 
(pred1p <- predict(mod1, type="response"))
exp(pl)/(1+exp(pl))
with(infarto,plot(prop , pred1p, asp=1)); abline(0,1)
 
(pred1y <- with(infarto, (S+N)*pred1p)) 
cbind(infarto$S, pred1y)
plot(infarto$S, pred1y, asp=1); abline(0,1)
 
 
with(infarto, plot(prop/(1-prop) ~ CK))
with(infarto, lines(pred1 ~ CK, type="b"))
with(infarto, plot(prop/(1-prop), pred1)); abline(0,1)
 
with(infarto, plot(prop ~ CK))
with(infarto, lines(pred1y ~ CK, type="b", col=2, pch=19, cex=0.5))
with(infarto, plot(prop, pred1y)); abline(0,1)
 
mod2 <- glm(cbind(S,N) ~ CK + I(CK^2), family=binomial, data=infarto)
mod3 <- glm(cbind(S,N) ~ CK + I(CK^2) + I(CK^3), family=binomial, data=infarto)
mod12 <- glm(cbind(S,N) ~ as.factor(CK), family=binomial, data=infarto)
 
(pred12p <- predict(mod12, type="response"))
infarto$prop 
 
 
c(logLik(mod0), logLik(mod1), logLik(mod2), logLik(mod3), logLik(mod12))
anova(mod0, mod1, mod2, mod3, mod12)
anova(mod0, mod1, mod2, mod3, mod12, test="Chisq")
 
mod1b <- glm(cbind(S,N) ~ CK, family=binomial(link="probit"), data=infarto)
mod1c <- glm(cbind(S,N) ~ CK, family=binomial(link="cauchit"), data=infarto)
logLik(mod1a)
logLik(mod1b)
logLik(mod1c)
 
 
mod3a <- glm(cbind(S,N) ~ poly(CK, 3), family=binomial, data=infarto)
mod3
mod3a
logLik(mod3)
logLik(mod3a)
deviance(mod3)
deviance(mod3a)
summary(mod3)
summary(mod3a)
anova(mod3)
anova(mod3a)
 
mod3b <- glm(cbind(S,N) ~ poly(CK, 3), family=binomial(link="probit"), data=infarto)
mod3c <- glm(cbind(S,N) ~ poly(CK, 3), family=binomial(link="cauchit"), data=infarto)
 
logLik(mod3a)
logLik(mod3b)
logLik(mod3c)
 
deviance(mod0)
-2*(logLik(mod0) - logLik(mod12))
 
deviance(mod1)
-2*(logLik(mod1) - logLik(mod12))
 
deviance(mod2)
-2*(logLik(mod2) - logLik(mod12))
 
deviance(mod3)
-2*(logLik(mod3) - logLik(mod12))
 
require(mgcv)
modgam <- gam(cbind(S,N) ~ s(CK), family=binomial, data=infarto)
modgam
logLik(modgam)
deviance(modgam)
 
plot(modgam)
 
with(infarto, plot(predict(modgam) ~ CK))
with(infarto, points(log(prop/(1-prop)) ~ CK))
## Fazer graficos de ajuste para o modelo 3 !!!

–

Considere o conjunto de dados "gala" com pacote "faraway"

require(faraway)
data(gala)
head(gala)
help(gala)

Vamos utilizar a função glm() do R, e seguem algumas sugestões (não exaustivas!!!) de tópicos a serem explorados pelas análises.

1. Qual modelo (GLM) poderia ser adotado?
2. Que gráficos exploratórios poderiam ser feitos?
3. O modelo NULO é suficiente para explicar os dados?
4. O modelo COMPLETO superior ao NULO para explicar os dados?
5. Como podemos avaliar a qualidade do ajuste do modelo completo?
6. Como podemos selecionar as covariáveis realmente relevantes?
7. Como obter intervalos para os parâmetros?
8. Como obter, avaliar e fazer gráficos do modelo predito?
9. Como avaliar diferentes opções para função de ligação?
10. Como avaliar se a linearidade das covariáveis selecionadas? (por exemplo um modelo com log(Area) seria melhor?)

Exemplo 1
A tabela abaixo apresenta dados de um estudo sobre "Acreditar em Vida Após a Morte" retirados de Wood (2006). O interesse aqui é utilizar estes dados para investigar se a crença está associada com o sexo. Embora várias abordagens e testes estatísticos sejam possíveis, nosso interesse investigar o problema a partir da especificação de GLM. Procure propor um GLM para o estudo, interpretar o(s) modelo(s) proposto(s), o nulo, o saturado, interpretando as quantidades de interesse nas análises e ainda comparando com resultados de outro(s) procedimento(s).

             Acredita
        -------------------
Sexo      Sim        Não
------------------------------
F         435	     147
M         375	     134
------------------------------

M <- cbind(c(435, 375), c(147, 134))
dimnames(M) <- list(c("F", "M"), c("S", "N"))
M
 
chisq.test(M)
 
addmargins(M)
Mesp <- outer(rowSums(M), colSums(M))/sum(M)
Mesp
 
(Chi2 <- sum(((M - Mesp)^2)/Mesp))
chisq.test(M, correct=F)
 
 
vam <- data.frame(Y=as.vector(M), Sexo = rownames(M), Acredita=rep(colnames(M), each=2))
vam
 
## Modelo 0
 
## Modelo 1: (independência)
## E(Y) = mu = n * Sexo * Acredita 
## log(mu)  = log(n) + log(Sexo) + log(Acredita) 
mod1 <- glm(Y ~ Sexo + Acredita, family=poisson(link="log"), data=vam)
model.matrix(mod1)
mod1
 
fitted(mod1)
Mesp
 
resid(mod1, type="pearson")
(M - Mesp)/sqrt(Mesp)
 
sum(resid(mod1, type="pearson")^2)
 
 
resid(mod1, type="deviance")
sum(resid(mod1, type="deviance")^2)
 
 
## Modelo 2: não independência
mod2 <- glm(Y ~ Sexo * Acredita, family=poisson, data=vam)
model.matrix(mod2)
mod2
 
fitted(mod2)
M
 
anova(mod1, mod2, test="Chisq")
 
## extensível a várias dimensões - modelos log-lineares
 
# sob link canônico X'y = X'\hat{mu}
#  - marginais preditas iguais a observadas em mod1

Exemplo 2
Novos casos de AIDS na Bélgica de 1981 a 1993.
Modelo Epidêmico (Venables & Ripley, MASS).
Interesse: A taxa de novos casos está reduzindo?

aids <- data.frame(
        y = c(12,14,33,50,67,74,123,141,165,204,253,246,240),
        ano = 1:13
)
aids
with(aids, plot(y ~ I(1980 + ano)))

Considere um modelo em que o número esperado de novos casos por ano seja:

Ajuste o modelo fazendo suposições necessárias e adequadas.
Procure avaliar a qualidade de ajuste e possíveis formas de tentar melhorar o ajuste com os dados disponíveis.
Represente o(s) modelo(s) ajustados graficamente (na escala das observações) procurando incluir no gráfico a incerteza das previsões

# mu_i = a . exp(b*ano_i)
# log(mu_i) = log(a) +  b*ano_i  = \beta_0 +  \beta_1 * ano_i
 
m1 <- glm(y ~ ano, family=poisson, data=aids) 
m1
summary(m1)
 
plot(resid(m1) ~fitted(m1))
lines(lowess(resid(m1) ~fitted(m1)))
 
par(mfrow=c(2,2))
plot(m1)
par(mfrow=c(1,1))
 
## tentando melhorar o ajuste
m2 <- glm(y ~ ano + I(ano^2), family=poisson, data=aids) 
m2
summary(m2)
anova(m1,m2, test="Chisq")
 
# nao faz sentido aqui mas é util com variaveis diferentes 
drop1(m2)
drop1(m2, test="Chisq")
 
ano.seq <- seq(1,13, l=100)
ypred <- predict(m2, data.frame(ano=ano.seq), se=TRUE)
str(ypred)
with(aids, plot(y ~ I(1980 + ano), ylim=c(0,280)))
lines(exp(ypred$fit) ~ I(ano.seq+1980))
lines(exp(ypred$fit - 2*ypred$se.fit) ~ I(ano.seq+1980), lty=2)
lines(exp(ypred$fit + 2*ypred$se.fit) ~ I(ano.seq+1980), lty=2)

1. Obter as funções de verossimilhança para os quatro modelos discutidos em aula para tabelas de contingência
2. proceder as análises com cada um dos modelos usando algum ambiente computacional
3. Examinar e justificar a razão pela qual podem ser ajustados com GLM com família Poisson
4. Comparar e discutir semelhanças e diferenças entre o ajuste de GLM com outros procedimentos (teste exato de fisher, teste chi-quadrado)

1. Mostrar como quantidades de interesse (predição) podem ser obtidas. Utilizar os exemplos vistos em aula.
   1. Mostrar como obter as predições na escala da resposta e da função de ligação
   2. Mostrar como obter os intervalos de predição
   3. Verificar os resultados com os retornados pela função predict() do R
   4. Algumas sugestões:
      1. como calcular os valores da curva de valores ajustados nos exemplos da creatinina vs infarto e no de "novos casos de AIDS na Bélgica". Ainda neste exemplo, como estimar a dose associada a uma certa probabilidade fixada de morte?
      2. como calcular as contagens esperadas no exemplo de crença vc sexo ?

Comandos do exemplo discutido em aula

## carregando o conjunto de dados DHF99 do pacote epicalc
require(epicalc)
data(DHF99)
head(DHF99)
help(DHF99)
summary(DHF99)
summary(DHF99)
## a variável village está como numérica mas de fato é um fator
DHF99$village <- as.factor(DHF99$village)
summary(DHF99)
##
## Ajustando GLM's
glm1 <- glm(containers ~ village + education, data=DHF99, family=poisson)
glm1
summary(glm1)
anova(glm1)
##
## vamos usar agora apenas a informação do tipo de vila
glm2 <- glm(containers ~ viltype + education, data=DHF99, family=poisson)
glm2
summary(glm1)
anova(glm1)
##
## comparando os ajustes
anova(glm2, glm1, test="Chisq")
## Portanto a informação individual de cada vila é relevante.
## Entretanto, como exemplo vamos supor que não dispomos da informação individual
## e apenas o tipo de vila.
## Neste caso o ajuste não é bom e vamos tentar alternativas
## 1. Modelo com interação
glm3 <- glm(containers ~ viltype * education, data=DHF99, family=poisson)
anova(glm2, glm3, test="Chisq")
anova(glm2, glm3, test="F")
## 2. Modelo Binomial negativo
require(MASS)
glm2BN <- glm.nb(containers ~ viltype + education, data=DHF99)
glm2BN
c(poisson=logLik(glm2), BN=logLik(glm2BN))
c(poisson=deviance(glm2), BN=deviance(glm2BN))
c(poisson=AIC(glm2), BN=AIC(glm2BN))
## comparando os coeficientes e erros padrão
summary(glm2)
summary(glm2BN)
## 3. Modelo quasipoisson
glm2Q <-
 glm(containers ~ viltype + education, data=DHF99, family=quasipoisson)
## comparando os coeficientes e erros padrão
summary(glm2)
summary(glm2Q)
## Note que os ajustes poderiam ser avaliados em mais detalhes
par(mfrow=c(2,2))
plot(glm1)
plot(glm2)
plot(glm3)
plot(glm2BN)
##
## As análises poderiam prosseguir de diversas formas. Alguns exemplos:
## - avaliando ainda outros modelos como por exemplo inflacionados de de zeros, hurdle .
## - verificando a relevância das covariáveis (e.g education)
## - Caso as vilas fossem consideradas poderia-se avaliar modelos mistos como vilas como efeito aleatórios