Modelos de Regressão para Dados de Contagem com o R

github.com/leg-ufpr/MRDCr

Dados referentes a um experimento realizado com o objetivo de investigar o efeito de uma intervenção, por parte do cuidador, no comportamento de ovelhas.

Para isso, foram consideradas ovelhas de duas linhagens distintas (pouco ou muito reativas), submetidas a dois tipos diferentes de intervenção (observação ou observação + intervenção).

A variável resposta aqui considerada é o número de mudanças na postura corporal do animal ao longo do período de observação (3 minutos).

Verificação do conteúdo e a estrutura dos dados.

str(datapost4)
## 'data.frame':    38 obs. of  3 variables:
##  $ Nposturas : num  3 11 15 6 1 1 4 16 2 1 ...
##  $ tratamento: Factor w/ 2 levels "Observ","Observ + Interv": 1 1 1 1 2 2 2 2 2 2 ...
##  $ linhagem  : Factor w/ 2 levels "Pouco reativo",..: 2 2 2 2 1 1 1 1 1 1 ...
summary(datapost4)
##    Nposturas                tratamento          linhagem 
##  Min.   : 0.000   Observ         :18   Pouco reativo:21  
##  1st Qu.: 3.000   Observ + Interv:20   Muito reativo:17  
##  Median : 6.000                                          
##  Mean   : 7.895                                          
##  3rd Qu.:11.750                                          
##  Max.   :35.000

Análise descritiva

tab <- data.frame(with(datapost4, table(tratamento, Nposturas)))

myColours <- brewer.pal(2,"Blues")
## Warning in brewer.pal(2, "Blues"): minimal value for n is 3, returning requested palette with 3 different levels
my.settings <- list(
  superpose.polygon=list(col=myColours[2:5], border="transparent"),
  strip.background=list(col=myColours[6]),
  strip.border=list(col="black")
)

bwplot(Nposturas ~ linhagem | tratamento,
         data=datapost4, 
         main="Mudanças de postura vs tratamento e linhagem", 
         xlab="Linhagem", ylab="Frequência")

### Variância e média amostrais por tratamento e linhagem.

mdp <- aggregate(Nposturas ~ tratamento:linhagem, datapost4, function(x) c(mean = mean(x), var = var(x)))
mdp
##        tratamento      linhagem Nposturas.mean Nposturas.var
## 1          Observ Pouco reativo       9.600000     28.266667
## 2 Observ + Interv Pouco reativo       4.636364     20.454545
## 3          Observ Muito reativo      12.125000    101.267857
## 4 Observ + Interv Muito reativo       6.222222     38.444444

Regressão poisson com estimação por máxima verossimilhança.

ajusteps <- glm(Nposturas ~ tratamento + linhagem, data=datapost4, family=poisson)
summary(ajusteps)
## 
## Call:
## glm(formula = Nposturas ~ tratamento + linhagem, family = poisson, 
##     data = datapost4)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.3579  -1.5592  -0.4320   0.7746   5.2884  
## 
## Coefficients:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                2.25083    0.09286  24.238  < 2e-16 ***
## tratamentoObserv + Interv -0.69665    0.12053  -5.780 7.48e-09 ***
## linhagemMuito reativo      0.25513    0.11549   2.209   0.0272 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 198.24  on 37  degrees of freedom
## Residual deviance: 158.46  on 35  degrees of freedom
## AIC: 298.65
## 
## Number of Fisher Scoring iterations: 5
### Avaliando possível efeito de interação.
ajustepsint <- glm(Nposturas ~ tratamento * linhagem, data=datapost4, family=poisson)
anova(ajusteps,ajustepsint, test = 'Chisq')
## Analysis of Deviance Table
## 
## Model 1: Nposturas ~ tratamento + linhagem
## Model 2: Nposturas ~ tratamento * linhagem
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1        35     158.46                     
## 2        34     158.40  1 0.063318   0.8013
exp(coef(ajusteps)[2])
## tratamentoObserv + Interv 
##                 0.4982513
### Estima-se que a taxa de variação de postura no grupo com intervenção seja aproximadamente
### a metade em relação ao grupo sem intervenção.

##### Estimação do parâmetro de dispersão. 
X2 <- sum(resid(ajusteps,type='pearson')**2)
phichap <- X2/ajusteps$df.residual
phichap
## [1] 5.094182

Diagnóstico do ajuste (gráficos).

##### Diagnóstico do modelo - gráficos padrão do R.
par(mfrow=c(2,2))
plot(ajusteps, pch = 20, cex = 1.25)

##### Gráfico quantil-quantil com envelopes simulados.
par(mfrow=c(1,1))
envelope(ajusteps)

Ajustando modelos por quase-verossimilhança.

### Modelo quasi poisson (V(mu) = mu).
ajuste12 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = 'quasipoisson')
# summary(ajuste12)

### Forma alternativa de declarar o Modelo quase-poisson (V(mu) = mu).
ajuste13 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = quasi(variance = 'mu', link='log'))
# summary(ajuste13)

### Modelo de quase-verossimilhança (V(mu) = mu^2).
ajuste14 <- glm(r4 ~ tratamento+linhagem, data=datapost4, family = quasi(variance = 'mu^2', link='log'))
summary(ajuste14)
## 
## Call:
## glm(formula = r4 ~ tratamento + linhagem, family = quasi(variance = "mu^2", 
##     link = "log"), data = datapost4)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5797  -0.6427  -0.1600   0.2356   1.6205  
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 2.2477     0.2337   9.616 2.34e-11 ***
## tratamentoObserv + Interv  -0.7007     0.2744  -2.554   0.0152 *  
## linhagemMuito reativo       0.2655     0.2756   0.964   0.3419    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasi family taken to be 0.7133735)
## 
##     Null deviance: 27.88  on 37  degrees of freedom
## Residual deviance: 22.67  on 35  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 5
### Gráficos de diagnóstico para o modelo de quase-verossimilhança.
par(mfrow = c(2,2))
plot(ajuste14, pch = 20, cex = 1.25)

### Gráficos quantil-quantil para os resíduos dos modelos Poisson e de Quase-Verossimilhança.
par(mfrow=c(1,2))
qqnorm(resid(ajusteps,type='deviance'), pch = 20, main = 'Poisson', las = 1)
qqline(resid(ajusteps,type='deviance'))
qqnorm(resid(ajuste14,type='deviance'), pch = 20, main = 'QL', las = 1)
qqline(resid(ajuste14,type='deviance'))

Usando estimação robusta e bootstrap.

estrb <- coeftest(ajusteps, vcov=sandwich) ### Estimador sanduíche
estrb ### Estimação robusta.
## 
## z test of coefficients:
## 
##                           Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)                2.25083    0.16483 13.6557 < 2.2e-16 ***
## tratamentoObserv + Interv -0.69665    0.26594 -2.6196  0.008804 ** 
## linhagemMuito reativo      0.25513    0.25452  1.0024  0.316144    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### Usando bootstrap (R=999 simulações)
ajusteboot <- Boot(ajusteps)
summary(ajusteboot) ### Resultados obtidos via bootstrap.
##                             R original   bootBias  bootSE  bootMed
## (Intercept)               999  2.25083 -0.0246026 0.18153  2.24489
## tratamentoObserv + Interv 999 -0.69665  0.0163349 0.27435 -0.69031
## linhagemMuito reativo     999  0.25513 -0.0098342 0.27313  0.24548

Resumo geral dos resultados.

erroz <- rbind(summary(ajusteps)$coefficients[2,2:3], summary(ajuste13)$coefficients[2,2:3], 
               summary(ajuste14)$coefficients[2,2:3], estrb[2,2:3], c(summary(ajusteboot)[2,4],
               mean(ajusteboot$t[,2]/summary(ajusteboot)[2,4])))

ics <- rbind(confint.default(ajusteps)[2,],confint.default(ajuste13)[2,], confint.default(ajuste14)[2,], 
estrb[2,1] + c(-1.96,1.96) * estrb[2,2], mean(ajusteboot$t[,2])+c(-1.96,1.96)*summary(ajusteboot)[2,4]) 

quadres <- cbind(erroz, ics)

rownames(quadres) <- c('Poisson', 'Quasi(mu)', 'Quasi(mu^2)', 'Robusto (sanduiche)', 'Bootstrap')

### Quadro resumo para as estimativas produzidas pelos quatro modelos para o efeito de intervenção.
kable(quadres, format = "markdown", caption = "Comparativo dos modelos ajustados")
Std. Error z value 2.5 % 97.5 %
Poisson 0.1205303 -5.779881 -0.9328858 -0.4604157
Quasi(mu) 0.2720405 -2.560835 -1.2298404 -0.1634612
Quasi(mu^2) 0.2744137 -2.553533 -1.2385656 -0.1628836
Robusto (sanduiche) 0.2659421 -2.619558 -1.2178973 -0.1754043
Bootstrap 0.2743507 -2.479731 -1.2180432 -0.1425886

Verificando efeito das observações com maiores resíduos na análise.

dadosexclud <- datapost4[-c(8,18,28),]

### Ajustando o modelo Poisson sem as três observações.
ajusteexcludpoi <- glm(Nposturas ~ tratamento+linhagem, data=dadosexclud, family = poisson)

### Estimativa do parâmetro de dispersão.
sum(resid(ajusteexcludpoi, type = 'pearson')**2)/ajusteexcludpoi$df.residual 
## [1] 1.990952
### Gráficos de diagnóstico.
par(mfrow=c(2,2))
plot(ajusteexcludpoi)

### Agora, o modelo quase-poisson
ajusteexcludquasi <- glm(Nposturas ~ tratamento+linhagem, data=dadosexclud, family = quasi(variance = 'mu', link='log'))

### Estimativas produzidas pelo modelo quasipoisson com e sem as três observações.

c1 <- compareCoefs(ajusteps, ajuste13, ajusteexcludpoi, ajusteexcludquasi, print = FALSE)
colnames(c1) <- c('Est. Poisson c/out','E.P. Poisson c/out',
                  'Est. Quasi c/out','E.P. Quasi c/out',
                   'Est. Poisson s/out','E.P. Poisson s/out',
                  'Est. Quasi s/out','E.P. Quasi s/out') 
kable(round(c1,3), align = 'c')
Est. Poisson c/out E.P. Poisson c/out Est. Quasi c/out E.P. Quasi c/out Est. Poisson s/out E.P. Poisson s/out Est. Quasi s/out E.P. Quasi s/out
(Intercept) 2.251 0.093 2.251 0.210 2.227 0.097 2.227 0.137
tratamentoObserv + Interv -0.697 0.121 -0.697 0.272 -0.886 0.144 -0.886 0.204
linhagemMuito reativo 0.255 0.115 0.255 0.261 0.005 0.134 0.005 0.190 
### Efeito da intervenção desconsiderando as três observações.
exp(coef(ajusteexcludpoi)[2])
## tratamentoObserv + Interv 
##                 0.4123804
### O efeito de intervenção aumenta, e torna-se mais significativo, mediante exclusão dos outliers.

Modelos de Regressão para Dados de Contagem com o R

61 RBRAS Departamento de Estatística - UFPR
23 a 25 de Maio de 2016 LEG, PET Estatística
Salvador - Bahia github.com/leg-ufpr/MRDCr