|
Modelos de Regressão para Dados de Contagem com o R
|
Walmes M. Zeviani, Eduardo E. Ribeiro Jr & Cesar A. Taconeli
library(MRDCr)
Implentando a função de log-verossimilhança do modelo COM-Poisson, definida como:
\[ \sum_i^n y_i \log(\lambda_i) - \nu \sum_i^n \log(y_i!) - \sum_i^n \log(Z(\lambda_i, \nu)) \]
em que \(Z(\lambda_i, \nu) = \sum_{j=0}^\infty \lambda_i^j (j!)^{-\nu}\) e \(\lambda_i = \exp(X_i\beta)\)
llcmp
## function (params, y, X, offset = NULL, sumto = ceiling(max(y)^1.5))
## {
## betas <- params[-1]
## phi <- params[1]
## nu <- exp(phi)
## if (is.null(offset))
## offset <- 0
## Xb <- X %*% betas + offset
## i <- 0:sumto
## zs <- sapply(Xb, function(loglambda) sum(exp(i * loglambda -
## nu * lfactorial(i))))
## Z <- sum(log(zs))
## ll <- sum(y * Xb - nu * lfactorial(y)) - Z
## return(-ll)
## }
## <environment: namespace:MRDCr>
Detalhes computacionais
Reparametrização do parâmetro \(\nu\) para \(\phi = \exp(\nu)\). Assim o espaçõ paramétrico do modelo são os reais \(\Re^n\).
Inclusão do termo offset. Somado diretamente ao preditor \(X_i \beta\), pois \(X_i \beta\) representa o parâmetro \(\lambda\) de locação, da distribuição COM-Poisson.
Truncamento da série infinita \(Z(\lambda_i)\). sumto
é tomado como argumento da função, que por padrão assume \(\max(\underline{y})^{1.2}\).
Para o cálculo de \(Z(\lambda_i)\) faz-se, minimizando problemas de overflow \[ \sum_{j=0}^\infty \lambda_i^j (j!)^{-\nu} = \sum_{j=0}^\infty \exp \left ( \log \left( \lambda_i^j (j!)^{-\nu} \right ) \right ) = \sum_{j=0}^\infty \exp(i \log(\lambda_i) - \nu \log(i!)) \]
Framework implementado em R que utiliza a forma de escrita de preditores no estilo de fórmulas, similar as funções lm
, glm
.
cmp
## function (formula, data, start = NULL, sumto = NULL, ...)
## {
## frame <- model.frame(formula, data)
## terms <- attr(frame, "terms")
## y <- model.response(frame)
## X <- model.matrix(terms, frame)
## off <- model.offset(frame)
## if (is.null(sumto))
## sumto <- ceiling(max(y)^1.5)
## if (is.null(start)) {
## m0 <- glm.fit(x = X, y = y, offset = off, family = poisson())
## start <- c(phi = 0, m0$coefficients)
## }
## bbmle::parnames(llcmp) <- names(start)
## model <- bbmle::mle2(llcmp, start = start, data = list(y = y,
## X = X, offset = off, sumto = sumto), vecpar = TRUE, ...)
## return(model)
## }
## <environment: namespace:MRDCr>
Um exemplo de como são construídas as matrizes, definidos os chutes iniciais e ajustados os modelos na função:
set.seed(2016)
x <- rep(1:3, 2)
t <- rnorm(6, 5)
y <- rpois(6, lambda = x*t)
(da <- data.frame(x, t, y))
## x t y
## 1 1 4.085258 2
## 2 2 6.001248 16
## 3 3 4.943577 13
## 4 1 5.296645 3
## 5 2 2.208529 5
## 6 3 4.717260 9
## Definindo o prditor do modelo
formula <- y ~ x + I(x^2) + offset(log(t))
##-------------------------------------------
## O framework
## Constrói as matrizes para ajuste do modelo
frame <- model.frame(formula, data = da)
(X <- model.matrix(formula, data = da))
## (Intercept) x I(x^2)
## 1 1 1 1
## 2 1 2 4
## 3 1 3 9
## 4 1 1 1
## 5 1 2 4
## 6 1 3 9
## attr(,"assign")
## [1] 0 1 2
(y <- model.response(frame))
## 1 2 3 4 5 6
## 2 16 13 3 5 9
(o <- model.offset(frame))
## [1] 1.4073849 1.7919674 1.5980892 1.6670736 0.7923267 1.5512280
## Utiliza como valores iniciais as estimativas dos parametros de um
## modelo GLM-Poisson
m0 <- glm.fit(x = X, y = y, offset = o, family = poisson())
start <- c(phi = 0, m0$coefficients)
## Otimiza a função de log-verossimilhança via bbmle
library(bbmle)
## Loading required package: stats4
parnames(llcmp) <- names(start)
mle2(llcmp, start = start,
data = list(y = y, X = X, offset = o, sumto = 50),
vecpar = TRUE)
##
## Call:
## mle2(minuslogl = llcmp, start = start, data = list(y = y, X = X,
## offset = o, sumto = 50), vecpar = TRUE)
##
## Coefficients:
## phi (Intercept) x I(x^2)
## 0.3338908 -4.5188932 5.4556688 -1.1173870
##
## Log-likelihood: -11.13
data(capdesfo)
str(capdesfo)
## 'data.frame': 125 obs. of 4 variables:
## $ est : Factor w/ 5 levels "vegetativo","botão floral",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ des : num 0 0 0 0 0 0.25 0.25 0.25 0.25 0.25 ...
## $ rept: int 1 2 3 4 5 1 2 3 4 5 ...
## $ ncap: int 10 9 8 8 10 11 9 10 10 10 ...
## help(capdesfo)
Experimento conduzido sob delineamento inteiramente casualizado em casa de vegetação onde avaliou-se o número de capulhos produzidos por plantas de algodão submetidas à níveis de desfolha artificial de remoção foliar em combinação com o estágio fenológico no qual a desfolha foi aplicada.
est
: Estágio fenológico com cinco níveis (vegetativo, botão floral, florecimento, maça, capulho);des
: Nível de desfolha artificial de remoção foliar (0, 25, 50, 75, 100%);ncap
: Número de capulhos de algodão produzidos ao final da ciclo cultura.## Experimento balanceado
xtabs(~est + des, data = capdesfo)
## des
## est 0 0.25 0.5 0.75 1
## vegetativo 5 5 5 5 5
## botão floral 5 5 5 5 5
## florecimento 5 5 5 5 5
## maça 5 5 5 5 5
## capulho 5 5 5 5 5
library(lattice)
library(latticeExtra)
(xy <- xyplot(ncap ~ des | est,
data = capdesfo,
xlab = "Nível de desfolha artificial",
ylab = "Número de capulhos produzidos",
type = c("p", "g", "smooth"),
panel = panel.beeswarm,
r = 0.05))
## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(ncap ~ est + des, data = capdesfo,
FUN = function(x) c(mean = mean(x), var = var(x))))
## est des ncap.mean ncap.var
## 1 vegetativo 0.00 9.0 1.0
## 2 botão floral 0.00 8.4 1.3
## 3 florecimento 0.00 9.4 2.8
## 4 maça 0.00 8.6 1.8
## 5 capulho 0.00 10.0 4.0
## 6 vegetativo 0.25 10.0 0.5
## 7 botão floral 0.25 9.4 3.3
## 8 florecimento 0.25 5.8 1.2
## 9 maça 0.25 7.0 2.0
## 10 capulho 0.25 10.0 3.5
## 11 vegetativo 0.50 8.6 0.8
## 12 botão floral 0.50 8.8 0.7
## 13 florecimento 0.50 5.8 1.7
## 14 maça 0.50 8.8 2.2
## 15 capulho 0.50 8.2 1.2
## 16 vegetativo 0.75 8.0 1.0
## 17 botão floral 0.75 8.8 2.7
## 18 florecimento 0.75 6.0 2.5
## 19 maça 0.75 6.2 0.2
## 20 capulho 0.75 8.8 1.2
## 21 vegetativo 1.00 6.2 1.2
## 22 botão floral 1.00 7.2 0.2
## 23 florecimento 1.00 4.6 0.3
## 24 maça 1.00 3.0 0.5
## 25 capulho 1.00 9.0 1.0
xlim <- ylim <- extendrange(c(mv$ncap), f = 0.05)
xyplot(ncap[, "var"] ~ ncap[, "mean"],
data = mv,
xlim = xlim,
ylim = ylim,
ylab = "Variância Amostral",
xlab = "Média Amostral",
panel = function(x, y) {
panel.xyplot(x, y, type = c("p", "r"), grid = TRUE)
panel.abline(a = 0, b = 1, lty = 2)
})
## Preditores considerados
f1 <- ncap ~ 1
f2 <- ncap ~ des + I(des^2)
f3 <- ncap ~ est:des + I(des^2)
f4 <- ncap ~ est:(des + I(des^2))
## Ajustando os modelos Poisson
m1P <- glm(f1, data = capdesfo, family = poisson)
m2P <- glm(f2, data = capdesfo, family = poisson)
m3P <- glm(f3, data = capdesfo, family = poisson)
m4P <- glm(f4, data = capdesfo, family = poisson)
## Ajustando os modelos COM-Poisson
m1C <- cmp(f1, data = capdesfo)
m2C <- cmp(f2, data = capdesfo)
m3C <- cmp(f3, data = capdesfo)
m4C <- cmp(f4, data = capdesfo)
## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P, m3P, m4P), logLik),
"COM-Poisson" = sapply(list(m1C, m2C, m3C, m4C), logLik))
## Poisson COM-Poisson
## [1,] -279.9328 -272.4794
## [2,] -271.3543 -256.0893
## [3,] -258.6741 -220.1977
## [4,] -255.8031 -208.2500
## Teste de razão de verossimilhanças
anova(m1P, m2P, m3P, m4P, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: ncap ~ 1
## Model 2: ncap ~ des + I(des^2)
## Model 3: ncap ~ est:des + I(des^2)
## Model 4: ncap ~ est:(des + I(des^2))
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 124 75.514
## 2 122 58.357 2 17.157 0.0001881 ***
## 3 118 32.997 4 25.360 4.258e-05 ***
## 4 114 27.255 4 5.742 0.2192611
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1C, m2C, m3C, m4C)
## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)
## Model 2: m2C, [llcmp]: phi+(Intercept)+des+I(des^2)
## Model 3: m3C, [llcmp]: phi+(Intercept)+I(des^2)+
## estvegetativo:des+estbotão floral:des+
## estflorecimento:des+estmaça:des+estcapulho:des
## Model 4: m4C, [llcmp]: phi+(Intercept)+estvegetativo:des+
## estbotão floral:des+estflorecimento:des+estmaça:des+
## estcapulho:des+estvegetativo:I(des^2)+estbotão
## floral:I(des^2)+estflorecimento:I(des^2)+
## estmaça:I(des^2)+estcapulho:I(des^2)
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 2 544.96
## 2 4 512.18 32.780 2 7.619e-08 ***
## 3 8 440.40 71.783 4 9.537e-15 ***
## 4 12 416.50 23.895 4 8.382e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Estimativas dos parâmetros
summary(m4P)
##
## Call:
## glm(formula = f4, family = poisson, data = capdesfo)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.06976 -0.31729 -0.02347 0.34779 1.27069
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.189560 0.063333 34.572 <2e-16 ***
## estvegetativo:des 0.436859 0.515560 0.847 0.3968
## estbotão floral:des 0.289715 0.507751 0.571 0.5683
## estflorecimento:des -1.242481 0.603706 -2.058 0.0396 *
## estmaça:des 0.364871 0.565803 0.645 0.5190
## estcapulho:des 0.008949 0.503760 0.018 0.9858
## estvegetativo:I(des^2) -0.805215 0.583915 -1.379 0.1679
## estbotão floral:I(des^2) -0.487850 0.566394 -0.861 0.3891
## estflorecimento:I(des^2) 0.672837 0.680195 0.989 0.3226
## estmaça:I(des^2) -1.310346 0.672780 -1.948 0.0515 .
## estcapulho:I(des^2) -0.019970 0.552505 -0.036 0.9712
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 75.514 on 124 degrees of freedom
## Residual deviance: 27.255 on 114 degrees of freedom
## AIC: 533.61
##
## Number of Fisher Scoring iterations: 4
summary(m4C)
## Maximum likelihood estimation
##
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y,
## X = X, offset = off, sumto = sumto), vecpar = TRUE)
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## phi 1.584744 0.127631 12.4166 < 2.2e-16 ***
## (Intercept) 10.896662 1.404311 7.7594 8.531e-15 ***
## estvegetativo:des 2.018744 1.140780 1.7696 0.0767909 .
## estbotão floral:des 1.343071 1.109168 1.2109 0.2259408
## estflorecimento:des -5.750451 1.479861 -3.8858 0.0001020 ***
## estmaça:des 1.594963 1.229254 1.2975 0.1944575
## estcapulho:des 0.037655 1.088012 0.0346 0.9723918
## estvegetativo:I(des^2) -3.724525 1.341954 -2.7754 0.0055125 **
## estbotão floral:I(des^2) -2.264721 1.254650 -1.8051 0.0710650 .
## estflorecimento:I(des^2) 3.134664 1.504387 2.0837 0.0371891 *
## estmaça:I(des^2) -5.894333 1.611937 -3.6567 0.0002555 ***
## estcapulho:I(des^2) -0.090110 1.193331 -0.0755 0.9398080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 416.4999
## Um dos problemas computacionais do modelo COM-Poisson é a obtenção da
## constante de normalização Z. Assim uma visualização pós ajuste para
## verificar se o ajuste proporcionou uma densidade válida se faz
## necessária
convergencez(m4C)
## Dado que o modelo COM-Poisson leva as mesmas estimativas pontuais que
## o modelo Poisson a análise de diagnóstico padrão pode ser utilizada
par(mfrow = c(2, 2))
plot(m4P)
##-------------------------------------------
## Testando a nulidade do parâmetro phi
## Usando o ajuste Poisson
trv <- 2 * (logLik(m4C) - logLik(m4P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
## [1] 95.10632 0.00000
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m4Cfixed <- cmp(f4, data = capdesfo, fixed = list("phi" = 0))
anova(m4C, m4Cfixed)
## Likelihood Ratio Tests
## Model 1: m4C, [llcmp]: phi+(Intercept)+estvegetativo:des+
## estbotão floral:des+estflorecimento:des+estmaça:des+
## estcapulho:des+estvegetativo:I(des^2)+estbotão
## floral:I(des^2)+estflorecimento:I(des^2)+
## estmaça:I(des^2)+estcapulho:I(des^2)
## Model 2: m4Cfixed, [llcmp]: (Intercept)+estvegetativo:des+
## estbotão floral:des+estflorecimento:des+estmaça:des+
## estcapulho:des+estvegetativo:I(des^2)+estbotão
## floral:I(des^2)+estflorecimento:I(des^2)+
## estmaça:I(des^2)+estcapulho:I(des^2)
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 12 416.50
## 2 11 511.61 95.106 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Via perfil de log-verossimilhança
perf <- profile(m4C, which = 1)
confint(perf)
## 2.5 % 97.5 %
## 1.323326 1.824516
plot(perf)
##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m4C)
Corr <- cov2cor(Vcov)
library(corrplot)
corrplot.mixed(Corr, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
## Predição pontual
pred <- with(capdesfo,
expand.grid(
est = levels(est),
des = seq(min(des), max(des), l = 20)
))
##-------------------------------------------
## Considerando a Poisson
mediasP <- exp(predict(m4P, newdata = pred))
aux <- data.frame(modelo = "Poisson", fit = mediasP)
predP <- cbind(pred, aux)
##-------------------------------------------
## Considerando a COM-Poisson
f4; f4[[2]] <- NULL; f4
## ncap ~ est:(des + I(des^2))
## ~est:(des + I(des^2))
X <- model.matrix(f4, data = pred)
## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m4C)[-1]
phi <- coef(m4C)[1]
loglambdas <- X %*% betas
## Aplicando a "inversa da função de ligação", ou seja, obtendo as
## contagens médias
mediasC <- sapply(loglambdas, FUN = function(p) {
y <- 0:50; py <- dcmp(y, p, phi, sumto = 100)
sum(y*py)
})
aux <- data.frame(modelo = "COM-Poisson", fit = mediasC)
predC <- cbind(pred, aux)
##-------------------------------------------
## Visualizando os valores preditos pelos dois modelos
da <- rbind(predP, predC)
update(xy, type = c("p", "g")) +
as.layer(xyplot(fit ~ des | est,
groups = modelo,
data = da, type = "l"))
## Predição intervalar
qn <- qnorm(0.975) * c(lwr = -1, upr = 1)
##-------------------------------------------
## Considerando a Poisson
aux <- predict(m4P, newdata = pred, se.fit = TRUE)
aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*")))
predP <- cbind(predP, aux)
##-------------------------------------------
## Considerando a COM-Poisson
## Obtendo os erros padrão das estimativas
## Obs.: Deve-se usar a matriz de variâncias e covariâncias
## condicional, pois os parâmetros de locação (betas) e dispersão
## (phi) não são ortogonais.
Vc <- Vcov[-1, -1] - Vcov[-1, 1] %*% solve(Vcov[1, 1]) %*% Vcov[1, -1]
U <- chol(Vc)
se <- sqrt(apply(X %*% t(U), MARGIN = 1, FUN = function(x) {
sum(x^2)
}))
aux <- c(loglambdas) + outer(se, qn, FUN = "*")
aux <- apply(aux, MARGIN = 2,
FUN = function(col) {
sapply(col, FUN = function(p) {
y <- 0:50; py <- dcmp(y, p, phi, sumto = 100)
sum(y*py)
})
})
predC <- cbind(predC, aux)
##-------------------------------------------
## Visualizando os valores preditos intervalares pelos dois modelos
da <- rbind(predP, predC)
update(xy, type = c("p", "g")) +
as.layer(xyplot(fit ~ des | est,
groups = modelo,
data = da,
type = "l",
ly = da$lwr,
uy = da$upr,
cty = "bands",
alpha = 0.3,
prepanel = prepanel.cbH,
panel.groups = panel.cbH,
panel = panel.superpose))
data(capmosca)
str(capmosca)
## 'data.frame': 60 obs. of 8 variables:
## $ dexp : int 0 0 0 0 0 0 0 0 0 0 ...
## $ vaso : Factor w/ 5 levels "1","2","3","4",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ planta : Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ alt : num 70.5 67.5 72 70 64 64.5 67.5 61 70 70 ...
## $ pesocap: num 25.2 NA 28.5 NA 31.8 ...
## $ nerep : int 4 4 3 5 5 5 5 5 4 5 ...
## $ ncap : int 4 4 2 5 5 5 5 5 4 5 ...
## $ nnos : int 13 14 14 15 15 11 14 11 14 15 ...
## help(capmosca)
Experimento conduzido sob delineamento inteiramente casualizado na Universidade Federal da Grande Dourados, cujo objetivo foi avaliar os impactos da exposição de plantas de algodão à alta infestação da praga mosca-branca. No experimento avaliou-se duas plantas por vaso, nesta análise tomaremos como unidade amostral o vaso e o interesse será somente na variável número de capulhos produzidos.
capmosca <- aggregate(ncap ~ vaso + dexp, data = capmosca, FUN = sum)
str(capmosca)
## 'data.frame': 30 obs. of 3 variables:
## $ vaso: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ dexp: int 0 0 0 0 0 1 1 1 1 1 ...
## $ ncap: int 8 7 10 10 9 5 7 8 8 11 ...
Assim as variáveis consideradas são definidas como:
dexp
: Dias de exposição à alta infestação de mosca-branca;ncap
: Número de capulhos de algodão produzidos ao final do experimento.## Experimento balanceado
xtabs(~dexp, data = capmosca)
## dexp
## 0 1 2 3 4 5
## 5 5 5 5 5 5
(xy <- xyplot(ncap ~ dexp,
data = capmosca,
xlab = "Dias de infestação",
ylab = "Número de capulhos produzidos",
type = c("p", "g", "smooth"),
panel = panel.beeswarm,
r = 0.05))
## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(ncap ~ dexp, data = capmosca,
FUN = function(x) c(mean = mean(x), var = var(x))))
## dexp ncap.mean ncap.var
## 1 0 8.8 1.7
## 2 1 7.8 4.7
## 3 2 6.8 5.2
## 4 3 6.8 1.7
## 5 4 7.4 2.3
## 6 5 7.6 2.3
## Preditores considerados
f1 <- ncap ~ 1
f2 <- ncap ~ dexp
f3 <- ncap ~ dexp + I(dexp^2)
## Ajustando os modelos Poisson
m1P <- glm(f1, data = capmosca, family = poisson)
m2P <- glm(f2, data = capmosca, family = poisson)
m3P <- glm(f3, data = capmosca, family = poisson)
## Ajustando os modelos COM-Poisson
m1C <- cmp(f1, data = capmosca)
m2C <- cmp(f2, data = capmosca)
m3C <- cmp(f3, data = capmosca)
## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P, m3P), logLik),
"COM-Poisson" = sapply(list(m1C, m2C, m3C), logLik))
## Poisson COM-Poisson
## [1,] -63.76981 -58.77209
## [2,] -63.52399 -58.13540
## [3,] -62.93611 -56.49274
## Teste de razão de verossimilhanças
anova(m1P, m2P, m3P, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: ncap ~ 1
## Model 2: ncap ~ dexp
## Model 3: ncap ~ dexp + I(dexp^2)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 29 11.993
## 2 28 11.501 1 0.49164 0.4832
## 3 27 10.325 1 1.17575 0.2782
anova(m1C, m2C, m3C)
## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)
## Model 2: m2C, [llcmp]: phi+(Intercept)+dexp
## Model 3: m3C, [llcmp]: phi+(Intercept)+dexp+I(dexp^2)
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 2 117.54
## 2 3 116.27 1.2734 1 0.2591
## 3 4 112.98 3.2853 1 0.0699 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Estimativas dos parâmetros
summary(m3P)
##
## Call:
## glm(formula = f3, family = poisson, data = capmosca)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.72211 -0.32208 0.01768 0.34877 1.13445
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.17571 0.13842 15.718 <2e-16 ***
## dexp -0.16923 0.13586 -1.246 0.213
## I(dexp^2) 0.02870 0.02638 1.088 0.277
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 11.993 on 29 degrees of freedom
## Residual deviance: 10.325 on 27 degrees of freedom
## AIC: 131.87
##
## Number of Fisher Scoring iterations: 4
summary(m3C)
## Maximum likelihood estimation
##
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y,
## X = X, offset = off, sumto = sumto), vecpar = TRUE)
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## phi 1.125858 0.261906 4.2987 1.718e-05 ***
## (Intercept) 6.824427 1.817389 3.7551 0.0001733 ***
## dexp -0.499138 0.266245 -1.8747 0.0608297 .
## I(dexp^2) 0.084627 0.050231 1.6848 0.0920312 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 112.9855
## Um dos problemas computacionais do modelo COM-Poisson é a obtenção da
## constante de normalização Z. Assim uma visualização pós ajuste para
## verificar se o ajuste proporcionou uma densidade válida se faz
## necessária
convergencez(m3C)
## Dado que o modelo COM-Poisson leva as mesmas estimativas pontuais que
## o modelo Poisson a análise de diagnóstico padrão pode ser utilizada
par(mfrow = c(2, 2))
plot(m3P)
##-------------------------------------------
## Testando a nulidade do parâmetro phi
## Usando o ajuste Poisson
trv <- 2 * (logLik(m3C) - logLik(m3P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
## [1] 12.88675 0.00033
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m3Cfixed <- cmp(f3, data = capmosca, fixed = list("phi" = 0))
anova(m3C, m3Cfixed)
## Likelihood Ratio Tests
## Model 1: m3C, [llcmp]: phi+(Intercept)+dexp+I(dexp^2)
## Model 2: m3Cfixed, [llcmp]: (Intercept)+dexp+I(dexp^2)
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 4 112.98
## 2 3 125.87 12.887 1 0.0003309 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Via perfil de log-verossimilhança
perf <- profile(m3C, which = 1)
confint(perf)
## 2.5 % 97.5 %
## 0.5622126 1.5977357
plot(perf)
##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m3C)
Corr <- cov2cor(Vcov)
library(corrplot)
corrplot.mixed(Corr, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
## Predição pontual/intervalar
pred <- with(capmosca,
expand.grid(
dexp = seq(min(dexp), max(dexp), l = 50)
))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
##-------------------------------------------
## Considerando a Poisson
aux <- predict(m3P, newdata = pred, se.fit = TRUE)
aux <- with(aux, exp(fit + outer(se.fit, qn, FUN = "*")))
aux <- data.frame(modelo = "Poisson", aux)
predP <- cbind(pred, aux)
##-------------------------------------------
## Considerando a COM-Poisson
f3; f3[[2]] <- NULL; f3
## ncap ~ dexp + I(dexp^2)
## ~dexp + I(dexp^2)
X <- model.matrix(f3, data = pred)
## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m3C)[-1]
phi <- coef(m3C)[1]
loglambdas <- X %*% betas
## Obtendo os erros padrão das estimativas
## Obs.: Deve-se usar a matriz de variâncias e covariâncias
## condicional, pois os parâmetros de locação (betas) e dispersão
## (phi) não são ortogonais.
Vc <- Vcov[-1, -1] - Vcov[-1, 1] %*% solve(Vcov[1, 1]) %*% Vcov[1, -1]
U <- chol(Vc)
se <- sqrt(apply(X %*% t(U), MARGIN = 1, FUN = function(x) {
sum(x^2)
}))
aux <- c(loglambdas) + outer(se, qn, FUN = "*")
aux <- apply(aux, MARGIN = 2,
FUN = function(col) {
sapply(col, FUN = function(p) {
y <- 0:30; py <- dcmp(y, p, phi, sumto = 50)
sum(y*py)
})
})
aux <- data.frame(modelo = "COM-Poisson", aux)
predC <- cbind(pred, aux)
##-------------------------------------------
## Visualizando os valores preditos intervalares pelos dois modelos
da <- rbind(predP, predC)
update(xy, type = c("p", "g")) +
as.layer(xyplot(fit ~ dexp,
groups = modelo,
data = da,
type = "l",
ly = da$lwr,
uy = da$upr,
cty = "bands",
alpha = 0.3,
prepanel = prepanel.cbH,
panel.groups = panel.cbH,
panel = panel.superpose))
data(ninfas)
str(ninfas)
## 'data.frame': 240 obs. of 8 variables:
## $ data : Date, format: "2009-12-11" ...
## $ dias : int 0 0 0 0 0 0 0 0 0 0 ...
## $ cult : Factor w/ 10 levels "BRS 239","BRS 243 RR",..: 3 3 3 3 2 2 2 2 4 4 ...
## $ bloco: Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ nsup : int 168 30 28 6 31 4 3 48 100 54 ...
## $ nmed : int 71 46 36 1 22 6 12 76 50 15 ...
## $ ninf : int 39 22 7 31 28 10 27 13 34 7 ...
## $ ntot : int 278 98 71 38 81 20 42 137 184 76 ...
## help(ninfas)
Experimento conduzido em casa de vegetação sob o delineamento de blocos casualizados. No experimento foram avaliadas plantas de diferentes cultivares de soja contabilizando o número de ninfas de mosca-branca nos folíolos dos terços superior, médio e inferior das plantas. As avaliações ocorreram em 6 datas dentre os 38 dias do estudo.
Nesta análise serão consideradas somente as cultivares com prefixo , sendo o número total de ninfas de mosca-branca nos folíolos a variável de interesse.
## Somente as cultivares que contém BRS na identificação
ninfas <- droplevels(subset(ninfas, grepl("BRS", x = cult)))
## Categorizando a variável dias em aval
ninfas$aval <- factor(ninfas$dias)
str(ninfas)
## 'data.frame': 96 obs. of 9 variables:
## $ data : Date, format: "2009-12-11" ...
## $ dias : int 0 0 0 0 0 0 0 0 0 0 ...
## $ cult : Factor w/ 4 levels "BRS 239","BRS 243 RR",..: 3 3 3 3 2 2 2 2 4 4 ...
## $ bloco: Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ nsup : int 168 30 28 6 31 4 3 48 100 54 ...
## $ nmed : int 71 46 36 1 22 6 12 76 50 15 ...
## $ ninf : int 39 22 7 31 28 10 27 13 34 7 ...
## $ ntot : int 278 98 71 38 81 20 42 137 184 76 ...
## $ aval : Factor w/ 6 levels "0","8","13","22",..: 1 1 1 1 1 1 1 1 1 1 ...
Assim as variáveis consideradas são definidas como:
## Experimento balanceado
xtabs(~aval + cult, data = ninfas)
## cult
## aval BRS 239 BRS 243 RR BRS 245 RR BRS 246 RR
## 0 4 4 4 4
## 8 4 4 4 4
## 13 4 4 4 4
## 22 4 4 4 4
## 31 4 4 4 4
## 38 4 4 4 4
(xy <- xyplot(ntot ~ aval | cult,
data = ninfas,
type = c("p", "g", "smooth"),
jitter.x = TRUE))
## Avaliando preliminarmente suposição de equidispersão
(mv <- aggregate(ntot ~ data + cult, data = ninfas,
FUN = function(x) c(mean = mean(x), var = var(x))))
## data cult ntot.mean ntot.var
## 1 2009-12-11 BRS 239 111.50000 475.66667
## 2 2009-12-19 BRS 239 149.00000 5262.66667
## 3 2009-12-24 BRS 239 170.25000 2889.58333
## 4 2010-01-02 BRS 239 54.75000 1634.91667
## 5 2010-01-11 BRS 239 6.25000 28.91667
## 6 2010-01-18 BRS 239 4.75000 40.91667
## 7 2009-12-11 BRS 243 RR 70.00000 2631.33333
## 8 2009-12-19 BRS 243 RR 76.25000 8435.58333
## 9 2009-12-24 BRS 243 RR 62.50000 2367.00000
## 10 2010-01-02 BRS 243 RR 30.00000 680.66667
## 11 2010-01-11 BRS 243 RR 4.50000 9.00000
## 12 2010-01-18 BRS 243 RR 3.50000 40.33333
## 13 2009-12-11 BRS 245 RR 121.25000 11522.25000
## 14 2009-12-19 BRS 245 RR 127.00000 9622.00000
## 15 2009-12-24 BRS 245 RR 113.50000 7052.33333
## 16 2010-01-02 BRS 245 RR 32.00000 568.66667
## 17 2010-01-11 BRS 245 RR 9.50000 51.00000
## 18 2010-01-18 BRS 245 RR 9.00000 134.66667
## 19 2009-12-11 BRS 246 RR 86.75000 4691.58333
## 20 2009-12-19 BRS 246 RR 124.00000 4790.66667
## 21 2009-12-24 BRS 246 RR 133.25000 5726.91667
## 22 2010-01-02 BRS 246 RR 38.75000 648.91667
## 23 2010-01-11 BRS 246 RR 7.00000 72.66667
## 24 2010-01-18 BRS 246 RR 4.50000 69.66667
xlim <- ylim <- extendrange(c(mv$ntot), f = 0.05)
xyplot(ntot[, "var"] ~ ntot[, "mean"],
data = mv,
xlim = xlim,
ylim = ylim,
ylab = "Variância Amostral",
xlab = "Média Amostral",
panel = function(x, y) {
panel.xyplot(x, y, type = c("p", "r"), grid = TRUE)
panel.abline(a = 0, b = 1, lty = 2)
})
## Preditores considerados
f1 <- ntot ~ bloco + cult + aval
f2 <- ntot ~ bloco + cult * aval
## Ajustando os modelos Poisson
m1P <- glm(f1, data = ninfas, family = poisson)
m2P <- glm(f2, data = ninfas, family = poisson)
## Ajustando os modelos COM-Poisson
m1C <- cmp(f1, data = ninfas, sumto = 600)
m2C <- cmp(f2, data = ninfas, sumto = 600)
## Warning in bbmle::mle2(llcmp, start = start, data = list(y = y, X
## = X, offset = off, : convergence failure: code=1 (iteration limit
## 'maxit' reached)
## Verossimilhancas dos modelos ajustados
cbind("Poisson" = sapply(list(m1P, m2P), logLik),
"COM-Poisson" = sapply(list(m1C, m2C), logLik))
## Poisson COM-Poisson
## [1,] -922.9782 -410.4440
## [2,] -879.2287 -407.1331
## Teste de razão de verossimilhanças
anova(m1P, m2P, test = "Chisq")
## Analysis of Deviance Table
##
## Model 1: ntot ~ bloco + cult + aval
## Model 2: ntot ~ bloco + cult * aval
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 84 1371.3
## 2 69 1283.8 15 87.499 2.897e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1C, m2C)
## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)+blocoII+blocoIII+
## blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
## aval8+aval13+aval22+aval31+aval38
## Model 2: m2C, [llcmp]: phi+(Intercept)+blocoII+blocoIII+
## blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
## aval8+aval13+aval22+aval31+aval38+cultBRS 243
## RR:aval8+cultBRS 245 RR:aval8+cultBRS 246 RR:aval8+
## cultBRS 243 RR:aval13+cultBRS 245 RR:aval13+cultBRS 246
## RR:aval13+cultBRS 243 RR:aval22+cultBRS 245 RR:aval22+
## cultBRS 246 RR:aval22+cultBRS 243 RR:aval31+cultBRS 245
## RR:aval31+cultBRS 246 RR:aval31+cultBRS 243 RR:aval38+
## cultBRS 245 RR:aval38+cultBRS 246 RR:aval38
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 13 820.89
## 2 28 814.27 6.6219 15 0.9673
## Estimativas dos parâmetros
summary(m1P)
##
## Call:
## glm(formula = f1, family = poisson, data = ninfas)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -13.2911 -3.0119 -0.7386 1.7743 12.3560
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.33397 0.03480 153.262 < 2e-16 ***
## blocoII -0.52019 0.03228 -16.114 < 2e-16 ***
## blocoIII -0.81268 0.03555 -22.857 < 2e-16 ***
## blocoIV -0.99360 0.03792 -26.204 < 2e-16 ***
## cultBRS 243 RR -0.69921 0.03894 -17.954 < 2e-16 ***
## cultBRS 245 RR -0.18595 0.03332 -5.582 2.38e-08 ***
## cultBRS 246 RR -0.23060 0.03373 -6.837 8.10e-12 ***
## aval8 0.20108 0.03416 5.887 3.94e-09 ***
## aval13 0.20788 0.03411 6.095 1.09e-09 ***
## aval22 -0.91822 0.04743 -19.360 < 2e-16 ***
## aval31 -2.65981 0.09908 -26.846 < 2e-16 ***
## aval38 -2.88525 0.11016 -26.191 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 7106.3 on 95 degrees of freedom
## Residual deviance: 1371.3 on 84 degrees of freedom
## AIC: 1870
##
## Number of Fisher Scoring iterations: 5
summary(m1C)
## Maximum likelihood estimation
##
## Call:
## bbmle::mle2(minuslogl = llcmp, start = start, data = list(y = y,
## X = X, offset = off, sumto = sumto), vecpar = TRUE)
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## phi -3.0829138 0.2092160 -14.7356 < 2.2e-16 ***
## (Intercept) 0.2434784 0.0518954 4.6917 2.709e-06 ***
## blocoII -0.0268210 0.0088711 -3.0234 0.0024994 **
## blocoIII -0.0428702 0.0113061 -3.7918 0.0001496 ***
## blocoIV -0.0533294 0.0129892 -4.1057 4.031e-05 ***
## cultBRS 243 RR -0.0380971 0.0113435 -3.3585 0.0007837 ***
## cultBRS 245 RR -0.0096912 0.0078249 -1.2385 0.2155260
## cultBRS 246 RR -0.0120559 0.0080379 -1.4999 0.1336482
## aval8 0.0103005 0.0079581 1.2943 0.1955470
## aval13 0.0106430 0.0079610 1.3369 0.1812566
## aval22 -0.0523118 0.0143159 -3.6541 0.0002581 ***
## aval31 -0.2305532 0.0452463 -5.0955 3.478e-07 ***
## aval38 -0.2708014 0.0530935 -5.1005 3.388e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 820.888
## Um dos problemas computacionais do modelo COM-Poisson é a obtenção da
## constante de normalização Z. Assim uma visualização pós ajuste para
## verificar se o ajuste proporcionou uma densidade válida se faz
## necessária
convergencez(m1C, incremento = 100, maxit = 10)
## Dado que o modelo COM-Poisson leva as mesmas estimativas pontuais que
## o modelo Poisson a análise de diagnóstico padrão pode ser utilizada
par(mfrow = c(2, 2))
plot(m1P)
##-------------------------------------------
## Testando a nulidade do parâmetro phi
## Usando o ajuste Poisson
trv <- 2 * (logLik(m1C) - logLik(m1P))
attributes(trv) <- NULL
round(c(trv, pchisq(trv, 1, lower = FALSE)), digits = 5)
## [1] 1025.068 0.000
## Reajustando o COM-Poisson para phi = 0 (ou equivalente nu = 1)
m1Cfixed <- cmp(f1, data = ninfas, fixed = list("phi" = 0))
anova(m1C, m1Cfixed)
## Likelihood Ratio Tests
## Model 1: m1C, [llcmp]: phi+(Intercept)+blocoII+blocoIII+
## blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
## aval8+aval13+aval22+aval31+aval38
## Model 2: m1Cfixed, [llcmp]: (Intercept)+blocoII+blocoIII+
## blocoIV+cultBRS 243 RR+cultBRS 245 RR+cultBRS 246 RR+
## aval8+aval13+aval22+aval31+aval38
## Tot Df Deviance Chisq Df Pr(>Chisq)
## 1 13 820.89
## 2 12 1845.96 1025.1 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Via perfil de log-verossimilhança
perf <- profile(m1C, which = 1)
confint(perf)
## 2.5 % 97.5 %
## -3.550614 -2.700415
plot(perf)
##-------------------------------------------
## Verificando a matriz ve variâncias e covariâncias
Vcov <- vcov(m1C)
Corr <- cov2cor(Vcov)
library(corrplot)
corrplot.mixed(Corr, lower = "number", upper = "ellipse",
diag = "l", tl.pos = "lt", tl.col = "black",
tl.cex = 0.8, col = brewer.pal(9, "Greys")[-(1:3)])
## Predição pontual/intervalar
pred <- with(ninfas,
expand.grid(
bloco = factor(levels(bloco)[1],
levels = levels(bloco)),
cult = levels(cult),
aval = levels(aval)
))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
f1; f1[[2]] <- NULL; f1
## ntot ~ bloco + cult + aval
## ~bloco + cult + aval
X <- model.matrix(f1, data = pred)
## Como não temos interesse na interpretação dos efeitos de blocos
## tomaremos a média desses efeitos para predição
bl <- attr(X, "assign") == 1
X[, bl] <- X[, bl] * 0 + 1/(sum(bl) + 1)
head(X)
## (Intercept) blocoII blocoIII blocoIV cultBRS 243 RR
## 1 1 0.25 0.25 0.25 0
## 2 1 0.25 0.25 0.25 1
## 3 1 0.25 0.25 0.25 0
## 4 1 0.25 0.25 0.25 0
## 5 1 0.25 0.25 0.25 0
## 6 1 0.25 0.25 0.25 1
## cultBRS 245 RR cultBRS 246 RR aval8 aval13 aval22 aval31 aval38
## 1 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0
## 3 1 0 0 0 0 0 0
## 4 0 1 0 0 0 0 0
## 5 0 0 1 0 0 0 0
## 6 0 0 1 0 0 0 0
library(multcomp)
##-------------------------------------------
## Considerando a Poisson
aux <- exp(confint(glht(m1P, linfct = X),
calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Poisson", aux)
predP <- cbind(pred, aux)
##-------------------------------------------
## Considerando a COM-Poisson
## Obtendo os parâmetros da distribuição (lambdas e phi)
betas <- coef(m1C)[-1]
phi <- coef(m1C)[1]
loglambdas <- X %*% betas
## Obtendo os erros padrão das estimativas
## Obs.: Deve-se usar a matriz de variâncias e covariâncias
## condicional, pois os parâmetros de locação (betas) e dispersão
## (phi) não são ortogonais.
Vc <- Vcov[-1, -1] - Vcov[-1, 1] %*% solve(Vcov[1, 1]) %*% Vcov[1, -1]
U <- chol(Vc)
se <- sqrt(apply(X %*% t(U), MARGIN = 1, FUN = function(x) {
sum(x^2)
}))
aux <- c(loglambdas) + outer(se, qn, FUN = "*")
aux <- apply(aux, MARGIN = 2,
FUN = function(col) {
sapply(col, FUN = function(p) {
y <- 0:400; py <- dcmp(y, p, phi, sumto = 600)
sum(y*py)
})
})
aux <- data.frame(modelo = "COM-Poisson", aux)
predC <- cbind(pred, aux)
##-------------------------------------------
## Visualizando os valores preditos intervalares pelos dois modelos
da <- rbind(predP, predC)
da <- da[order(da$cult, da$aval, da$modelo), ]
source(paste0("https://gitlab.c3sl.ufpr.br/leg/legTools/raw/",
"issue%2315/R/panel.segplot.by.R"))
update(xy, type = c("p", "g"), alpha = 0.5) +
as.layer(segplot(
aval ~ lwr + upr | cult, centers = fit,
data = da,
horizontal = FALSE, draw = FALSE, lwd = 2, grid = TRUE,
panel = panel.segplot.by, groups = modelo, f = 0.1,
pch = 1:nlevels(pred$modelo)+2, as.table = TRUE))
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 |