Reparametrization of COM-Poisson Regression Models with Applications in the Analysis of Experimental Count Data
Eduardo Elias Ribeiro Junior
Walmes Marques Zeviani
Wagner Hugo Bonat
Clarice Garcia Borges Demétrio
March 31, 2017
General functions
##======================================================================
## Tranformations
lambda2mu <- function(lambda, nu) {
phi <- log(nu)
mu <- lambda^(1/nu) - (nu - 1)/(2 * nu)
list("mu" = mu, "phi" = phi)
}
mu2lambda <- function(mu, phi) {
nu <- exp(phi)
lambda <- (mu + (nu-1)/(2*nu))^nu
list("lambda" = lambda, "nu" = nu)
}
##======================================================================
## Compute momentes mean and variance
calc_moments <- Vectorize(FUN = function(mu, phi) {
pars <- mu2lambda(mu, phi)
mu <- cmpreg::calc_mean_cmp(log(pars[["lambda"]]), phi)
va <- cmpreg::calc_var_cmp(log(pars[["lambda"]]), phi)
## mu <- do.call(compoisson::com.mean, pars)
## va <- do.call(compoisson::com.var, pars)
c("mean" = mu, "var" = va)
}, vectorize.args = c("mu", "phi"))
##======================================================================
## Negative log-likelihood functions
## COM-Poisson Model (Original parameterization)
llcmp <- function (params, y, X, offset = NULL, sumto = NULL) {
if (!is.null(offset)) {
stop("This model doesn't support offset yet.")
}
phi <- params[1]
nu <- exp(phi)
Xb <- X %*% params[-1]
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)
}
## COM-Poisson Model (mean-parameterization)
llcmp2 <- function (params, y, X, offset = NULL, sumto = NULL) {
if (!is.null(offset)) {
stop("This model doesn't support offset yet.")
}
phi2 <- params[1]
nu <- exp(phi2)
##-------------------------------------------
mu <- exp(X %*% params[-1])
Xb <- nu * log(mu + (nu - 1)/(2 * nu))
##-------------------------------------------
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)
}
##======================================================================
## Framework for fit count models (using bbmle::mle2 function)
fitcm <- function(formula, data, start = NULL,
model = c("CP", "CP2"), ...) {
dots <- list(...)
model <- match.arg(model)
##-------------------------------------------
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(start)) {
m0 <- glm.fit(x = X, y = y, offset = off, family = poisson())
start <- c(0, m0$coefficients)
}
##-------------------------------------------
if (model == "CP") {
names(start)[1] <- "phi"
bbmle::parnames(llcmp) <- names(start)
sumto <- dots$sumto
dots$sumto <- NULL
args <- append(
list(minuslog = llcmp, start = start,
data = list(y = y, X = X, offset = off,
sumto = sumto, model = model),
vecpar = TRUE), dots)
fit <- suppressWarnings(
do.call(what = bbmle::mle2, args = args))
}
if (model == "CP2") {
names(start)[1] <- "phi2"
bbmle::parnames(llcmp2) <- names(start)
sumto <- dots$sumto
dots$sumto <- NULL
args <- append(
list(minuslog = llcmp2, start = start,
data = list(y = y, X = X, offset = off,
sumto = sumto, model = model),
vecpar = TRUE), dots)
fit <- suppressWarnings(
do.call(what = bbmle::mle2, args = args))
}
slot(fit, "call.orig") <- match.call()
return(fit)
}
##======================================================================
## Analysis of Deviance Table to glm and mle2 objetcs
getAnova <- function(..., print = TRUE) {
model.list <- list(...)
test <- "Chisq"
cls <- sapply(model.list, function(x) class(x)[1])
if (all(cls == "glm")) {
fam <- sapply(model.list, function(x) x$family$family)
if (all(fam == "quasipoisson")) {
test <- "F"
}
}
nps <- sapply(model.list, function(x) attr(logLik(x), "df"))
aic <- sapply(model.list, function(x) AIC(x))
if (test == "Chisq") {
lls <- sapply(model.list, function(x) logLik(x))
cst <- 2 * diff(lls)
pvs <- pchisq(cst, df = diff(nps), lower.tail = FALSE)
tab <- cbind(
"np" = nps,
"ll" = lls,
"AIC" = aic,
"2*dif" = c(NA, cst),
"dif np" = c(NA, diff(nps)),
"Pr(>Chisq)" = c(NA, pvs))
}
if (test == "F") {
## ver pág. 187 expresssão 7.10 em Modern Applied, Ripley
sigma <- sapply(model.list, function(x) summary(x)$dispersion)
df.sigma <- sapply(model.list, function(x) x$df.residual)
dev <- sapply(model.list, function(x) deviance(x))
dfs <- c(NA, -diff(dev))
dnp <- c(NA, diff(nps))
F <- c(NA, -diff(dev))/(dnp * sigma)
pvs <- pf(F, df1 = dnp, df2 = df.sigma, lower.tail = FALSE)
tab <- cbind(
"np" = nps,
"dev" = dev,
"AIC" = aic,
"F" = F,
"dif np" = dnp,
"Pr(>F)" = pvs)
}
rownames(tab) <- 1:length(model.list)
if (print) printCoefmat(tab, na.print = "", cs.ind = 1)
invisible(tab)
}
## Get the estimated coefficients of the parameters
getCoefs <- function(model, digits = 3, rownames = NULL) {
##-------------------------------------------
glue <- function(est, d = digits) {
est
## apply(est, 1, FUN = function(x)
## paste0(format(x[1], digits = d, nsmall = d), " (",
## format(x[2], digits = d, nsmall = d), ")"))
}
##-------------------------------------------
cls <- class(model)[1]
est <- switch(
cls,
"glm" = {
fam <- model$family$family
if (!fam %in% c("poisson", "quasipoisson")) {
stop("Família de distribuições não contemplada")
}
switch(
fam,
"poisson" = {
glue(rbind(c(NA, NA),
summary(model)$coef[, c(1, 3)]))
},
"quasipoisson" = {
glue(rbind(c(summary(model)$dispersion, NA),
summary(model)$coef[, c(1, 3)]))
}
)
},
"mle2" = {
glue(bbmle::summary(model)@coef[, c(1, 3)])
}
)
if (!is.null(rownames)) {
rownames(est) <- rownames
}
colnames(est) <- c("Est", "Est/EP")
return(est)
}
##======================================================================
## Predict Values
## Get the interval values of linear predictors (for Delta method with
## cholesky decomposition)
cholV_eta <- function(V, X, b, qn) {
eta <- c(X %*% b)
if (length(qn) > 1) {
U <- chol(V)
stderr <- sqrt(
apply(X %*% t(U), MARGIN = 1, FUN = function(x) sum(x^2))
)
eta <- sweep(x = outer(stderr, qn, FUN = "*"),
MARGIN = 1,
STATS = eta,
FUN = "+")
}
return(eta)
}
## Get the means (inverse link function)
calc_mean <- function(eta, dispersion, model = c("CP", "CP2"),
tol = 1e-5, offset = NULL, ...) {
if (is.null(offset)) {
offset <- 1L
}
model <- match.arg(model)
names(eta) <- NULL
names(dispersion) <- NULL
pars <- data.frame(eta = eta, dispersion = dispersion)
##-------------------------------------------
if (model == "CP") {
##-------------------------------------------
pars <- transform(pars, lambda = exp(eta),
nu = exp(dispersion))
sumto <- list(...)$sumto
##-------------------------------------------
approxmu <- with(pars, lambda^(1/nu) - (nu - 1)/(2 * nu))
sigma <- with(pars, (1/nu) * approxmu)
sigma <- ifelse(sigma < 0, 0, sigma)
ymax <- ceiling(max(approxmu + 5 * sqrt(sigma)))
etamax <- max(pars$eta)
dispersionmin <- min(pars$dispersion)
pmax <- exp(
-llcmp(params = c(dispersionmin, etamax), y = ymax, X = 1,
sumto = sumto))
while (pmax > tol) {
ymax <- ymax + 1
pmax <- exp(
-llcmp(params = c(dispersionmin, etamax),
y = ymax, X = 1, sumto = sumto))
}
y <- 1:ymax
estmean <- mapply(
eta = pars$eta,
nu = pars$nu,
MoreArgs = list(y = y),
FUN = function(eta, nu, y) {
i <- 0:sumto
zs <- sapply(eta, function(eta) {
sum(exp(i * eta - nu * lfactorial(i)))
})
py <- exp(y * eta - nu * lfactorial(y) - sum(log(zs)))
sum(y * py)
},
SIMPLIFY = TRUE)
}
if (model == "CP2") {
estmean <- exp(eta)
}
names(estmean) <- NULL
return(c(estmean))
}
## Calculate the confidence interval of predict values
predictcm <- function(object, newdata, offset = NULL, level = 0.95) {
##-------------------------------------------
if (missing(newdata)) {
X <- object@data$X
} else {
X <- newdata
}
##-------------------------------------------
qn <- -qnorm((1 - level[1])/2)
qn <- qn * c(fit = 0, lwr = -1, upr = 1)
#--------------------------------------------
V <- vcov(object)
Vc <- V[-1, -1] - V[-1, 1] %*% solve(V[1, 1]) %*% V[1, -1]
## Vc <- V[-1, -1]
eta <- cholV_eta(Vc, X, b = coef(object)[-1], qn = qn)
##-------------------------------------------
args <- list(dispersion = coef(object)[1],
model = object@data$model,
offset = offset)
if (object@data$model == "CP") {
args$sumto <- object@data$sumto
}
pred <- apply(as.matrix(eta), MARGIN = 2,
FUN = function(x) {
do.call(what = calc_mean,
args = append(list(eta = x), args))
})
colnames(pred) <- names(qn)
return(pred)
}
##======================================================================
## Trellis panel functions for lattice graphics (use wzRfun)
library(wzRfun)
body(panel.cbH)[[14]][[3]][[3]] <-
quote(panel.lines(xo, lyo, col = col.line, lty = list(...)$lty))
body(panel.cbH)[[14]][[3]][[4]] <-
quote(panel.lines(xo, uyo, col = col.line, lty = list(...)$lty))
## To profile log-likelihood
myprofile <- function(model, which) {
sel <- c("param", "z", "focal")
prof <- profile(model, which = which)
as.data.frame(prof)[, sel]
}
## To plot profile in lattice
xyprofile <- function(prof, conf = c(0.9, 0.95, 0.99),
namestrip = NULL, subset = 4,
scales.x = "free", ...) {
##-------------------------------------------
conf <- conf[order(conf, decreasing = TRUE)]
da <- subset(prof, abs(z) <= subset)
##-------------------------------------------
fl <- levels(da$param)
if (!is.null(namestrip)) {
fl <- namestrip
}
xyplot(abs(z) ~ focal | param,
data = da,
layout = c(NA, 1),
## ylab = expression(abs(z)~~(sqrt(~Delta~"deviance"))),
ylab = expression(sqrt(2~(L(hat(phi))-L(phi)))),
scales = list(x = scales.x),
type = c("l", "g"),
strip = strip.custom(
factor.levels = fl),
panel = function(x, y, subscripts, ...) {
conf <- c(0.9, 0.95, 0.99)
hl <- sqrt(qchisq(conf, 1))
##-------------------------------------------
mle <- x[y == 0]
xl <- x[x < mle]; yl <- y[x < mle]
xr <- x[x > mle]; yr <- y[x > mle]
##-------------------------------------------
funleft <- approxfun(x = yl, y = xl)
funright <- approxfun(x = yr, y = xr)
vxl <- funleft(hl)
vxr <- funright(hl)
vz <- c(hl, hl)
##-------------------------------------------
panel.xyplot(x, y, ...)
panel.arrows(c(vxl, vxr), 0, c(vxl, vxr), vz,
code = 1, length = 0.1, lty = 2,
col = "gray40")
panel.segments(vxl, vz, vxr, vz, lty = 2,
col = "gray40")
panel.abline(h = 0, v = mle, lty = 3)
panel.text(x = rep(mle, 2), y = vz + 0.1,
labels = paste(conf*100, "%"),
col = "gray20")
}, ...)
}Graphics configuration
## Report settings
library(knitr)
opts_chunk$set(
warning = FALSE,
message = FALSE,
fig.width = 9,
fig.height = 4,
dev.args = list(family = "Palatino"))## Load the lattice packages and customize outputs graphics
library(lattice)
library(latticeExtra)
library(grid)
## http://www.magesblog.com/2013/04/how-to-change-alpha-value-of-colours-in.html
add.alpha <- function(col, alpha = 1){
apply(sapply(col, col2rgb)/255, 2,
function(x) rgb(x[1], x[2], x[3], alpha = alpha))
}
## Define Colors
mycol <- c(1, "#377EB8", "#E41A1C", "#4DAF4A",
"#ff00ff", "#FF7F00", "#984EA3", "#FFFF33", "#808080")
myreg <- colorRampPalette(c(mycol[3], "gray90", mycol[2]))(100)
## Trellis graphical style.
ps <- list(
superpose.symbol = list(
col = mycol, pch = 1,
fill = add.alpha(mycol, alpha = 0.4)),
box.rectangle = list(col = 1, fill = c("gray70")),
box.umbrella = list(col = 1, lty = 1),
box.dot = list(pch = "|"),
dot.symbol = list(col = 1, pch = 19),
dot.line = list(col = "gray50", lty = 3),
plot.symbol = list(col = 1),
plot.line = list(col = 1),
plot.polygon = list(col = "gray95"),
superpose.line = list(col = mycol, lty = 1),
superpose.polygon = list(col = mycol),
strip.background = list(col = c("gray90", "gray70")),
regions = list(col = myreg),
par.sub.text = list(
font = 1, just = "left", cex = 0.9,
x = grid::unit(10, "mm"))
)
trellis.par.set(ps)
## Alternative settings
ac <- list(pad1 = 0.5, pad2 = 0.5, tck = 0.5)
ps2 <- list(
layout.widths = list(
left.padding = 0.25,
right.padding = -1,
ylab.axis.padding = 0),
layout.heights = list(
bottom.padding = 0.25,
top.padding = 0,
axis.xlab.padding = 0,
xlab.top = 0),
axis.components = list(
bottom = ac, top = ac,
left = ac, right = ac)
)Moments Approximation
## =====================================================================
## Analysis of Approximation moments
## Eduardo Junior
## edujrrib@gmail.com
## 2017-03-16
## =====================================================================## Load package and functions
source("config.R")
source("functions.R")
## Colors for legends
cols <- trellis.par.get("superpose.line")$col## Convergence of Z(lambda, nu) constant
computeZ <- function(lambda, nu, maxit = 1e4, tol = 1e-5) {
z <- vector("numeric", maxit)
j = 1
z[j] <- exp(j * log(lambda) - nu * lfactorial(j))
while (abs(z[j] - 0) > tol && j <= maxit) {
j = j + 1
z[j] <- exp(j * log(lambda) - nu * lfactorial(j))
}
if (lambda == 1 & nu == 0) z <- Inf
return(sum(z, na.rm = TRUE) + 1)
}
grid <- expand.grid(
lambda = c(0.5, 1, 5, 10, 30, 50),
nu = seq(0, 1, length.out = 11))
grid$z <- apply(grid, 1, function(par) {
computeZ(lambda = par[1], nu = par[2], maxit = 1e4)
})
xt <- xtabs(z ~ nu + lambda, data = grid)
xt## lambda
## nu 0.5 1 5 10
## 0 1.999992e+00 Inf Inf Inf
## 0.1 1.916269e+00 7.642745e+00 Inf Inf
## 0.2 1.856385e+00 5.253523e+00 3.167527e+273 Inf
## 0.3 1.810601e+00 4.319695e+00 1.602324e+29 2.541984e+282
## 0.4 1.774097e+00 3.803204e+00 4.710590e+10 1.326456e+56
## 0.5 1.744133e+00 3.469505e+00 1.339958e+06 3.668680e+22
## 0.6 1.719005e+00 3.233784e+00 2.050181e+04 4.993080e+12
## 0.7 1.697585e+00 3.057328e+00 2.372940e+03 3.685683e+08
## 0.8 1.679084e+00 2.919745e+00 6.488009e+02 2.698812e+06
## 0.9 1.662937e+00 2.809192e+00 2.740935e+02 1.466388e+05
## 1 1.648721e+00 2.718282e+00 1.484132e+02 2.202647e+04
## lambda
## nu 30 50
## 0 Inf Inf
## 0.1 Inf Inf
## 0.2 Inf Inf
## 0.3 Inf Inf
## 0.4 Inf Inf
## 0.5 3.319764e+196 Inf
## 0.6 1.730117e+76 4.627475e+177
## 0.7 4.929488e+39 6.933351e+81
## 0.8 5.086211e+24 3.425421e+46
## 0.9 1.801582e+17 2.191499e+30
## 1 1.068647e+13 5.184706e+21## Study the approximation## Mean and variance relationship
aux <- expand.grid(
mu = seq(2, 30, length.out = 50),
phi = seq(log(0.3), log(2.5), length.out = 50))
moments <- mapply(FUN = calc_moments,
mu = aux$mu,
phi = aux$phi,
SIMPLIFY = FALSE)
grid <- cbind(aux, t(do.call(cbind, moments)))
grid <- transform(grid, va = mu / exp(phi))
col <- brewer.pal(n = 8, name = "RdBu")
myreg <- colorRampPalette(colors = col)
xy1 <- xyplot(var ~ mean,
groups = phi,
data = grid,
type = "l",
lwd = 2,
axis = axis.grid,
xlab = expression(E(X)),
ylab = expression(V(X)),
sub = "(a)",
legend = list(
top = list(
fun = draw.colorkey,
args = list(
key = list(
space = "top",
col = myreg(length(unique(grid$phi))),
at = unique(grid$phi),
draw = FALSE)))),
par.settings = modifyList(
ps2, list(superpose.line = list(
col = myreg(length(unique(grid$phi))))
)
),
panel = function(x, y, ...) {
panel.xyplot(x, y, ...)
panel.curve(1*x, min(x), max(x), lty = 2)
},
page = function(...) {
grid.text(expression(log(nu)),
just = "bottom",
x = unit(0.10, "npc"),
y = unit(0.83, "npc"))
})## Errors in approximations for E[Y] and V[Y]
grid <- transform(grid,
emu = (mu - mean)^2,
eva = (va - var)^2)
myreg <- colorRampPalette(c("gray90", "gray20"))(100)
xy2 <- levelplot(emu ~ phi + mu, data = grid,
aspect = "fill",
col.regions = myreg,
xlab = expression(phi),
ylab = expression(mu),
sub = "(b)",
colorkey = list(space = "top"),
par.settings = ps2,
panel = function(x, y, z, ...) {
panel.levelplot(x, y, z, ...)
panel.curve(10 - ( exp(x) - 1)/(2 * exp(x)),
lty = 2)
panel.abline(v = 0, lty = 2)
})
xy3 <- levelplot(eva ~ phi + mu, data = grid,
aspect = "fill",
col.regions = myreg,
xlab = expression(phi),
ylab = expression(mu),
sub = "(c)",
colorkey = list(space = "top"),
par.settings = ps2,
panel = function(x, y, z, ...) {
panel.levelplot(x, y, z, ...)
panel.curve((10 - ( exp(x) - 1)/
(2 * exp(x)))/exp(x), lty = 2)
panel.abline(v = 0, lty = 2)
})print(xy1, split = c(1, 1, 3, 1), more = TRUE)
print(xy2, split = c(2, 1, 3, 1), more = TRUE)
print(xy3, split = c(3, 1, 3, 1), more = FALSE)
## COM-Poisson probability mass function (mean parametrization)
dcmp <- Vectorize(FUN = function(y, mu, phi, sumto = 100)
exp(-llcmp2(c(phi, log(mu)), y = y, X = 1, sumto = sumto)),
vectorize.args = c("y", "mu", "phi"))
grid <- expand.grid(mu = c(2, 8, 15), phi = log(c(0.5, 1, 2.5)))
y <- 0:30
py <- mapply(FUN = dcmp,
mu = grid$mu,
phi = grid$phi,
MoreArgs = list(y = y, sumto = 100),
SIMPLIFY = FALSE)
grid <- cbind(grid[rep(1:nrow(grid), each = length(y)), ],
y = y,
py = unlist(py))## COM-Poisson p.m.f. to different combination betwenn phi and mu
leg_phi <- parse(
text = paste("phi == \"",
formatC(unique(grid$phi), 1, format = "f"),
"\""))
barchart(py ~ y | factor(mu),
groups = factor(phi),
data = grid,
horizontal = FALSE,
layout = c(NA, 1),
as.table = TRUE,
axis = axis.grid,
origin = 0,
xlim = extendrange(y, f = 0.01),
border = "transparent",
scales = list(x = list(at = pretty(y))),
ylab = expression(P(Y==y)),
xlab = expression(y),
par.settings = list(
superpose.polygon = list(col = c(mycol[2:4])),
superpose.line = list(col = c(mycol[2:4]))
),
auto.key = list(
columns = 3,
rectangles = FALSE,
lines = TRUE,
text = leg_phi
),
strip = strip.custom(
strip.names = TRUE,
var.name = expression(mu == ""),
sep = ""))
Defoliation experiment
## =====================================================================
## Analysis of Artificial Defoliation in Cotton Plants
## Eduardo Junior
## edujrrib@gmail.com
## 2017-03-16
## =====================================================================## Load package and functions
library(bbmle)
library(multcomp)
library(plyr)
source("config.R")
source("functions.R")
## Colors for legends
cols <- trellis.par.get("superpose.line")$col## Load data
## data(cottonBolls, package = "cmpreg")
cottonBolls <- read.table("./data/cottonBolls.txt",
header = TRUE, sep = "\t")
str(cottonBolls)## 'data.frame': 125 obs. of 4 variables:
## $ est : Factor w/ 5 levels "botão floral",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ 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 ...## Exploratory analysis
## Scatter plot with a beewarm taste
xy1 <- xyplot(ncap ~ des | est,
data = cottonBolls,
layout = c(2, 3),
as.table = TRUE,
type = c("p", "smooth", "g"),
xlab = "Níveis de desfolha artificial",
ylab = "Número de capulhos produzidos",
spread = 0.05,
panel = panel.beeswarm)
## Sample variance vs sample mean (evidence in favor of the
## underdispersion).
mv <- aggregate(ncap ~ est + des, data = cottonBolls,
FUN = function(x) c(mean = mean(x), var = var(x)))
xlim <- ylim <- extendrange(c(mv$ncap), f = 0.05)
xy2 <- xyplot(ncap[, "var"] ~ ncap[, "mean"],
data = mv,
type = c("p", "r", "g"),
xlim = xlim,
ylim = ylim,
xlab = expression("Média"~"amostral"~(bar(y))),
ylab = expression("Variância"~"amostral"~(s^2)),
panel = function(...) {
panel.xyplot(...)
panel.abline(a = 0, b = 1, lty = 2)
})print(xy1, split = c(1, 1, 2, 1), more = TRUE)
print(xy2, split = c(2, 1, 2, 1), more = FALSE)
## Fit models
mnames <- c("PO", "C1", "C2", "QP")
## Predictor, following Zeviani et al. (2014)
form1 <- ncap ~ est:(des + I(des^2))
m1PO <- glm(form1, data = cottonBolls, family = poisson)
system.time(
m1C1 <- fitcm(form1, data = cottonBolls, model = "CP", sumto = 50)
)## user system elapsed
## 4.656 0.000 4.655system.time(
m1C2 <- fitcm(form1, data = cottonBolls, model = "CP2", sumto = 50)
)## user system elapsed
## 3.516 0.000 3.518m1QP <- glm(form1, data = cottonBolls, family = quasipoisson)
models.ncap <- list(m1PO, m1C1, m1C2, m1QP)
names(models.ncap) <- mnames
## Numbers of calls to loglik and numerical gradient
models.ncap$C1@details$counts## function gradient
## 86 39models.ncap$C2@details$counts## function gradient
## 63 19## LRT and profile extra parameter
## LRT between Poisson and COM-Poisson (test: phi == 0)
getAnova(m1PO, m1C2)## np ll AIC 2*dif dif np Pr(>Chisq)
## 1 11 -255.8 533.61
## 2 12 -208.4 440.80 94.811 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1profs.ncap <- lapply(list(c(m1C1, "phi"), c(m1C2, "phi2")),
function(x) myprofile(x[[1]], x[[2]]))
profs.ncap <- do.call("rbind", profs.ncap)
## Plot profile extra parameter
snames <- parse(text = c("'COM-Poisson'~(phi)",
"'COM-Poisson'[mu]~(phi)"))
xyprofile(profs.ncap, namestrip = snames,
ylim = extendrange(c(0, 3)),
xlab = "Dispersion/Precision parameter")
## Goodness of fit measures and estimate parameters
## GoF measures
measures.ncap <- sapply(models.ncap, function(x)
c("LogLik" = logLik(x), "AIC" = AIC(x), "BIC" = BIC(x)))
measures.ncap## PO C1 C2 QP
## LogLik -255.8031 -208.2500 -208.3977 NA
## AIC 533.6062 440.4999 440.7953 NA
## BIC 564.7177 474.4397 474.7351 NA## Get the estimates
est <- lapply(models.ncap, FUN = function(x) getCoefs(x))
est.ncap <- do.call(cbind, est)
est.ncap## Est Est/EP Est Est/EP
## NA NA 1.58474438 12.41658734
## (Intercept) 2.189560352 34.57242139 10.89666236 7.75943954
## estbotão floral:des 0.289715410 0.57058563 1.34307061 1.21088166
## estcapulho:des 0.008949298 0.01776502 0.03765465 0.03460873
## estflorescimento:des -1.242480824 -2.05808949 -5.75045075 -3.88580567
## estmaça:des 0.364871430 0.64487355 1.59496264 1.29750502
## estvegetativo:des 0.436859418 0.84734944 2.01874374 1.76961737
## estbotão floral:I(des^2) -0.487850140 -0.86132683 -2.26472100 -1.80506216
## estcapulho:I(des^2) -0.019970497 -0.03614535 -0.09010993 -0.07551137
## estflorescimento:I(des^2) 0.672836705 0.98918234 3.13466390 2.08368190
## estmaça:I(des^2) -1.310346157 -1.94765883 -5.89433299 -3.65667768
## estvegetativo:I(des^2) -0.805215332 -1.37899472 -3.72452543 -2.77544938
## Est Est/EP Est
## 1.581686716 12.39219235 0.241024719
## (Intercept) 2.190453506 74.63972962 2.189560352
## estbotão floral:des 0.287644027 1.22265905 0.289715410
## estcapulho:des 0.007556596 0.03238496 0.008949298
## estflorescimento:des -1.247157801 -4.42015869 -1.242480824
## estmaça:des 0.350034297 1.32799559 0.364871430
## estvegetativo:des 0.434997229 1.81938426 0.436859418
## estbotão floral:I(des^2) -0.485805354 -1.84993286 -0.487850140
## estcapulho:I(des^2) -0.018939962 -0.07399356 -0.019970497
## estflorescimento:I(des^2) 0.678814707 2.13491360 0.672836705
## estmaça:I(des^2) -1.287531326 -4.09513789 -1.310346157
## estvegetativo:I(des^2) -0.803316413 -2.96133854 -0.805215332
## Est/EP
## NA
## (Intercept) 70.42048395
## estbotão floral:des 1.16222453
## estcapulho:des 0.03618553
## estflorescimento:des -4.19211766
## estmaça:des 1.31354143
## estvegetativo:des 1.72596409
## estbotão floral:I(des^2) -1.75443459
## estcapulho:I(des^2) -0.07362438
## estflorescimento:I(des^2) 2.01486318
## estmaça:I(des^2) -3.96718170
## estvegetativo:I(des^2) -2.80887111## Prediction
## Data for prediction
pred <- with(cottonBolls,
expand.grid(
est = levels(est),
des = seq(min(des), max(des), l = 20)
))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
## Design matrix for prediction
X <- model.matrix(update(form1, NULL~.), pred)
## Considering Poisson
aux <- exp(confint(
glht(m1PO, linfct = X), calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Poisson", aux)
predPO.ncap <- cbind(pred, aux)
## Considering COM-Poisson
aux <- predictcm(m1C1, newdata = X)
aux <- data.frame(modelo = "COM-Poisson", aux)
predC1.ncap <- cbind(pred, aux)
## Considering COM-Poisson (mean parametrization)
aux <- predictcm(m1C2, newdata = X)
aux <- data.frame(modelo = "COM-Poisson2", aux)
predC2.ncap <- cbind(pred, aux)
## Considering Quasi-Poisson
aux <- exp(confint(
glht(m1QP, linfct = X), calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Quasi-Poisson", aux)
predQP.ncap <- cbind(pred, aux)
## Representing the confidence intervals
pred.ncap <- rbind(predPO.ncap, predC1.ncap, predC2.ncap, predQP.ncap)
## Legend
key <- list(columns = 2,
lines = list(col = cols[1:4], lty = rev(1:4)),
text = list(parse(
text = c("'Poisson'", "'COM-Poisson'",
"'COM-Poisson'[mu]", "'Quasi-Poisson'"))
))## Graph
xyplot(ncap ~ des | est,
data = cottonBolls,
layout = c(NA, 1),
as.table = TRUE,
type = c("p", "g"),
xlab = "Níveis de desfolha artificial",
ylab = "Número de capulhos produzidos",
spread = 0.05,
alpha = 0.6, key = key,
panel = panel.beeswarm) +
as.layer(
xyplot(fit ~ des | est,
auto.key = TRUE,
data = pred.ncap,
groups = modelo,
type = "l",
layout = c(NA, 1),
as.table = TRUE,
ly = pred.ncap$lwr, uy = pred.ncap$upr,
cty = "bands", fill = "gray80", alpha = 0.1,
panel = panel.superpose,
panel.groups = panel.cbH,
prepanel = cmpreg::prepanel.cbH,
lty = rev(1:4))
)
## Correlation between estimates
corr.ncap <- purrr::map_df(models.ncap[c("C1", "C2")],
function(x) cov2cor(vcov(x))[1, -1])
corr.ncap## # A tibble: 11 × 2
## C1 C2
## <dbl> <dbl>
## 1 0.995247654 5.217831e-04
## 2 0.152633164 -2.282951e-04
## 3 0.004269042 -1.753199e-04
## 4 -0.489517939 -6.773895e-04
## 5 0.161353485 -1.484546e-03
## 6 0.222881120 -1.647535e-04
## 7 -0.227589625 1.799200e-04
## 8 -0.009463862 1.167071e-04
## 9 0.262871798 6.332628e-04
## 10 -0.457756799 1.774117e-03
## 11 -0.349629408 9.967896e-05Soybean experiment
## =====================================================================
## Analysis of number of grains in Soyabean under Soil Moisture and
## Potassium Doses
## Eduardo Junior
## edujrrib@gmail.com
## 2017-03-16
## =====================================================================## Load package and functions
library(bbmle)
library(multcomp)
library(plyr)
source("config.R")
source("functions.R")
## Colors for legends
cols <- trellis.par.get("superpose.line")$col## Load data
## data(soyaBeans, package = "cmpreg")
soyaBeans <- read.table("./data/soyaBeans.txt",
header = TRUE, sep = "\t")
soyaBeans$umid <- as.factor(soyaBeans$umid)
soyaBeans <- soyaBeans[-74, ] ## Incorrect observation
soyaBeans <- transform(soyaBeans, K = K/100)## Exploratory analysis
## Scatter plot
xy1 <- xyplot(ngra ~ K | umid,
data = soyaBeans,
xlab = "Nível de adubação potássica",
ylab = "Número de grãos por parcela",
type = c("p", "g", "smooth"),
as.table = TRUE,
layout = c(2, 2),
strip = strip.custom(
strip.names = TRUE, var.name = "Umidade",
factor.levels = paste0(levels(soyaBeans$umid), "%")))
## Sample variance vs sample mean (evidence in favor of the
## overdispersion).
mv <- aggregate(ngra ~ K + umid, data = soyaBeans,
FUN = function(x) c(mean = mean(x), var = var(x)))
xlim <- ylim <- extendrange(c(mv$ngra), f = 0.05)
xy2 <- xyplot(ngra[, "var"] ~ ngra[, "mean"],
data = mv,
type = c("p", "r", "g"),
xlim = xlim,
ylim = ylim,
xlab = expression("Média"~"amostral"~(bar(y))),
ylab = expression("Variância"~"amostral"~(s^2)),
panel = function(...) {
panel.xyplot(...)
panel.abline(a = 0, b = 1, lty = 2)
})print(xy1, split = c(1, 1, 2, 1), more = TRUE)
print(xy2, split = c(2, 1, 2, 1), more = FALSE)
## Fit models
## Predictor
form2 <- ngra ~ bloc + umid * K + I(K^2)
m2PO <- glm(form2, data = soyaBeans, family = poisson)
system.time(
m2C1 <- fitcm(form2, data = soyaBeans, model = "CP", sumto = 700)
)## user system elapsed
## 27.976 0.004 27.985system.time(
m2C2 <- fitcm(form2, data = soyaBeans, model = "CP2", sumto = 700)
)## user system elapsed
## 14.044 0.004 14.050m2QP <- glm(form2, data = soyaBeans, family = quasipoisson)
models.ngra <- list(m2PO, m2C1, m2C2, m2QP)
names(models.ngra) <- mnames
## Numbers of calls to loglik and numerical gradient
models.ngra$C1@details$counts## function gradient
## 264 72models.ngra$C2@details$counts## function gradient
## 78 20## LRT and profile extra parameter
## LRT between Poisson and COM-Poisson (test: phi == 0)
getAnova(m2PO, m2C2)## np ll AIC 2*dif dif np Pr(>Chisq)
## 1 11 -340.08 702.16
## 2 12 -325.23 674.47 29.697 1 5.052e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1profs.ngra <- lapply(list(c(m2C1, "phi"), c(m2C2, "phi2")),
function(x) myprofile(x[[1]], x[[2]]))
profs.ngra <- do.call("rbind", profs.ngra)
## Plot profile extra parameter
snames <- parse(text = c("'COM-Poisson'~(phi)",
"'COM-Poisson'[mu]~(phi)"))
xyprofile(profs.ngra, namestrip = snames,
ylim = extendrange(c(0, 3)),
xlab = "Dispersion/Precision parameter",
scales.x = NULL)
## Goodness of fit measures and estimate parameters
## GoF measures
measures.ngra <- sapply(models.ngra, function(x)
c("LogLik" = logLik(x), "AIC" = AIC(x), "BIC" = BIC(x)))
measures.ngra## PO C1 C2 QP
## LogLik -340.0818 -325.2409 -325.2334 NA
## AIC 702.1635 674.4817 674.4668 NA
## BIC 727.5083 702.1305 702.1156 NA## Get the estimates
est <- lapply(models.ngra, FUN = function(x) getCoefs(x))
est.ngra <- do.call(cbind, est)
est.ngra## Est Est/EP Est Est/EP Est
## NA NA -0.778547541 -4.7208417 -0.78208463
## (Intercept) 4.86663406 144.2885947 2.232040988 6.0415142 4.86659270
## blocII -0.01939070 -0.7302129 -0.008929678 -0.4939046 -0.01940469
## blocIII -0.03662632 -1.3732566 -0.016867366 -0.9212309 -0.03663866
## blocIV -0.10559332 -3.8890255 -0.048633814 -2.4223326 -0.10555415
## blocV -0.09164670 -3.2997067 -0.042207142 -2.1020159 -0.09169581
## umid50 0.13202198 3.6471269 0.060853865 2.2948915 0.13201627
## umid62.5 0.12430785 3.4318953 0.057301026 2.1772441 0.12430805
## K 0.61599366 11.0138894 0.283863485 4.7291342 0.61612811
## I(K^2) -0.27594121 -10.2501073 -0.127151306 -4.5890201 -0.27600680
## umid50:K 0.14555067 4.2679607 0.066986558 2.6137768 0.14556134
## umid62.5:K 0.16481872 4.8294234 0.075857611 2.8842900 0.16481527
## Est/EP Est Est/EP
## -4.7371043 2.61509049 NA
## (Intercept) 97.7807529 4.86663406 89.2254286
## blocII -0.4950102 -0.01939070 -0.4515503
## blocIII -0.9305848 -0.03662632 -0.8491968
## blocIV -2.6337588 -0.10559332 -2.4049022
## blocV -2.2366064 -0.09164670 -2.0404783
## umid50 2.4715036 0.13202198 2.2553166
## umid62.5 2.3257542 0.12430785 2.1222213
## K 7.4639917 0.61599366 6.8107878
## I(K^2) -6.9458464 -0.27594121 -6.3384790
## umid50:K 2.8921765 0.14555067 2.6392289
## umid62.5:K 3.2722880 0.16481872 2.9864271## Prediction
## Data for prediction
pred <- with(soyaBeans,
expand.grid(
bloc = factor(levels(bloc)[1], levels = levels(bloc)),
umid = levels(umid),
K = seq(min(K), max(K), l = 20)
))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
## Design matrix for prediction
X <- model.matrix(update(form2, NULL ~ .), pred)
bl <- attr(X, "assign") == 1
X[, bl] <- X[, bl] + 1/(sum(bl) + 1)
## Considering Poisson
aux <- exp(confint(
glht(m2PO, linfct = X), calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Poisson", aux)
predPO.ngra <- cbind(pred, aux)
## Considering COM-Poisson
aux <- predictcm(m2C1, newdata = X)
aux <- data.frame(modelo = "COM-Poisson", aux)
predC1.ngra <- cbind(pred, aux)
## Considering COM-Poisson (mean parametrization)
aux <- predictcm(m2C2, newdata = X)
aux <- data.frame(modelo = "COM-Poisson2", aux)
predC2.ngra <- cbind(pred, aux)
## Considering Quasi-Poisson
aux <- exp(confint(
glht(m2QP, linfct = X), calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Quasi-Poisson", aux)
predQP.ngra <- cbind(pred, aux)
## Representing the confidence intervals
pred.ngra <- rbind(predPO.ngra, predC1.ngra, predC2.ngra, predQP.ngra)
## Legend
key <- list(columns = 2,
lines = list(col = cols[1:4], lty = rev(1:4)),
text = list(parse(
text = c("'Poisson'", "'COM-Poisson'",
"'COM-Poisson'[mu]", "'Quasi-Poisson'"))
)
)## Graph
update(xy1, layout = c(NA, 1), type = c("p", "g"),
alpha = 0.6, key = key) +
as.layer(
xyplot(fit ~ K | umid,
data = pred.ngra,
groups = modelo,
type = "l",
ly = pred.ngra$lwr, uy = pred.ngra$upr,
cty = "bands", fill = "gray80", alpha = 0.1,
panel = panel.superpose,
panel.groups = panel.cbH,
prepanel = cmpreg::prepanel.cbH,
lty = rev(1:4))
)
## Correlation between estimates
corr.ngra <- purrr::map_df(models.ngra[c("C1", "C2")],
function(x) cov2cor(vcov(x))[1, -1])
corr.ngra## # A tibble: 11 × 2
## C1 C2
## <dbl> <dbl>
## 1 0.99807743 2.815933e-04
## 2 -0.08098072 9.509210e-05
## 3 -0.15103753 8.359179e-05
## 4 -0.39706227 -3.047646e-04
## 5 -0.34459432 3.487660e-04
## 6 0.37561348 1.943150e-05
## 7 0.35632417 -2.164181e-05
## 8 0.77459471 -4.902517e-04
## 9 -0.75177928 4.852495e-04
## 10 0.42905820 -5.339187e-05
## 11 0.47342759 3.661366e-05Nitrofen experiment
## =====================================================================
## Analysis of number of grains in Soyabean under Soil Moisture and
## Potassium Doses
## Eduardo Junior
## edujrrib@gmail.com
## 2017-03-16
## =====================================================================## Load package and functions
library(bbmle)
library(multcomp)
library(plyr)
source("config.R")
source("functions.R")
## Colors for legends
cols <- trellis.par.get("superpose.line")$col## Load data
## data(nitrofen, package = "boot")
## data(Paula, package = "labestData")
## nitrofen <- PaulaEx4.6.20
nitrofen <- read.table("./data/nitrofen.txt",
header = TRUE, sep = "\t")
nitrofen <- transform(nitrofen, dose = dose/100)## Exploratory analysis
## Scatter plot
xy1 <- xyplot(novos ~ dose,
data = nitrofen,
xlab = "Nível de adubação potássica",
ylab = "Número de grãos por parcela",
type = c("p", "g", "smooth"),
spread = 0.05,
panel = panel.beeswarm)
## Sample variance vs sample mean (evidence in favor of the
## relationship with dispersion and covariate).
mv <- aggregate(novos ~ dose, data = nitrofen,
FUN = function(x) c(mean = mean(x), var = var(x)))
xlim <- ylim <- extendrange(c(mv$novos), f = 0.05)
xy2 <- xyplot(novos[, "var"] ~ novos[, "mean"],
data = mv,
type = c("p", "r", "g"),
xlim = xlim,
ylim = ylim,
xlab = expression("Média"~"amostral"~(bar(y))),
ylab = expression("Variância"~"amostral"~(s^2)),
panel = function(...) {
panel.xyplot(...)
panel.abline(a = 0, b = 1, lty = 2)
})print(xy1, split = c(1, 1, 2, 1), more = TRUE)
print(xy2, split = c(2, 1, 2, 1), more = FALSE)
## Fit models
mnames <- c("PO", "C1", "C2", "QP")
## Predictors
form31 <- novos ~ dose
form32 <- novos ~ dose + I(dose^2)
form33 <- novos ~ dose + I(dose^2) + I(dose^3)
predictors <- list("pred1" = form31, "pred2" = form32, "pred3" = form33)
fmodels.ovos <- lapply(predictors, function(form) {
PO <- glm(form, data = nitrofen, family = poisson)
C1 <- fitcm(form, data = nitrofen, model = "CP", sumto = 100)
C2 <- fitcm(form, data = nitrofen, model = "CP2", sumto = 100)
QP <- glm(form, data = nitrofen, family = quasipoisson)
list("PO" = PO, "C1" = C1, "C2" = C2, "QP" = QP)
})## LRT for nested models
## Poisson
auxPO <- lapply(fmodels.ovos, function(x) x$PO)
do.call("getAnova", auxPO)## np ll AIC 2*dif dif np Pr(>Chisq)
## 1 2 -180.67 365.33
## 2 3 -147.01 300.02 67.3192 1 2.309e-16 ***
## 3 4 -144.09 296.18 5.8354 1 0.01571 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## COM-Poisson standard
auxC1 <- lapply(fmodels.ovos, function(x) x$C1)
do.call("getAnova", auxC1)## np ll AIC 2*dif dif np Pr(>Chisq)
## 1 3 -167.95 341.91
## 2 4 -146.96 301.93 41.9797 1 9.222e-11 ***
## 3 5 -144.06 298.13 5.7998 1 0.01603 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## COM-Poisson mean-parameterized
auxC2 <- lapply(fmodels.ovos, function(x) x$C2)
do.call("getAnova", auxC2)## np ll AIC 2*dif dif np Pr(>Chisq)
## 1 3 -167.65 341.30
## 2 4 -146.95 301.90 41.4047 1 1.238e-10 ***
## 3 5 -144.06 298.13 5.7725 1 0.01628 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Quasi-Poisson
auxQP <- lapply(fmodels.ovos, function(x) x$QP)
do.call("getAnova", auxQP)## np dev AIC F dif np Pr(>F)
## 1 2 123.929
## 2 3 56.610 60.840 1 5.073e-10 ***
## 3 4 50.774 5.659 1 0.02157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Separe the choose models
form3 <- form33
m3PO <- fmodels.ovos$pred3$PO
m3C1 <- fmodels.ovos$pred3$C1
m3C2 <- fmodels.ovos$pred3$C2
m3QP <- fmodels.ovos$pred3$QP
models.ovos <- list(m3PO, m3C1, m3C2, m3QP)
names(models.ovos) <- mnames
## Numbers of calls to loglik and numerical gradient
models.ovos$C1@details$counts## function gradient
## 56 10models.ovos$C2@details$counts## function gradient
## 46 8## LRT and profile extra parameter
## LRT between Poisson and COM-Poisson (test: phi == 0)
getAnova(m3PO, m3C2)## np ll AIC 2*dif dif np Pr(>Chisq)
## 1 4 -144.09 296.18
## 2 5 -144.06 298.13 0.052935 1 0.818profs.novos <- lapply(list(c(m3C1, "phi"), c(m3C2, "phi2")),
function(x) myprofile(x[[1]], x[[2]]))
profs.novos <- do.call("rbind", profs.novos)
## Plot profile extra parameter
snames <- parse(text = c("'COM-Poisson'~(phi)",
"'COM-Poisson'[mu]~(phi)"))
xyprofile(profs.novos, namestrip = snames,
ylim = extendrange(c(0, 3)),
xlab = "Dispersion/Precision parameter",
scales.x = NULL) +
layer(panel.abline(v = 0, col = 2, lty = 2))
## Goodness of fit measures and estimate parameters
## GoF measures
measures.ovos <- sapply(models.ovos, function(x)
c("LogLik" = logLik(x), "AIC" = AIC(x), "BIC" = BIC(x)))
measures.ovos## PO C1 C2 QP
## LogLik -144.0901 -144.0644 -144.0636 NA
## AIC 296.1802 298.1288 298.1272 NA
## BIC 303.8283 307.6889 307.6874 NA## Get the estimates
est <- lapply(models.ovos, FUN = function(x) getCoefs(x))
est.ovos <- do.call(cbind, est)
est.ovos## Est Est/EP Est Est/EP Est
## NA NA 0.04822702 0.2355371 0.04742685
## (Intercept) 3.47673043 62.8166664 3.64941764 4.8498534 3.47687492
## dose -0.08604628 -0.4327654 -0.09137570 -0.4475093 -0.08786760
## I(dose^2) 0.15290176 0.8633605 0.16121072 0.8782597 0.15470220
## I(dose^3) -0.09723110 -2.3977670 -0.10208018 -2.2293683 -0.09764724
## Est/EP Est Est/EP
## 0.2317296 1.03117215 NA
## (Intercept) 64.3083270 3.47673043 61.8599120
## dose -0.4522797 -0.08604628 -0.4261740
## I(dose^2) 0.8938402 0.15290176 0.8502107
## I(dose^3) -2.4635375 -0.09723110 -2.3612469## Prediction
## Data for prediction
pred <- with(nitrofen,
data.frame("dose" = seq(min(dose), max(dose),
length.out = 100)))
qn <- qnorm(0.975) * c(fit = 0, lwr = -1, upr = 1)
## Design matrix for prediction
X <- model.matrix(update(form3, NULL ~ .), pred)
## Considering Poisson
aux <- exp(confint(
glht(m3PO, linfct = X), calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Poisson", aux)
predPO.novos <- cbind(pred, aux)
## Considering COM-Poisson
aux <- predictcm(m3C1, newdata = X)
aux <- data.frame(modelo = "COM-Poisson", aux)
predC1.novos <- cbind(pred, aux)
## Considering COM-Poisson (mean parametrization)
aux <- predictcm(m3C2, newdata = X)
aux <- data.frame(modelo = "COM-Poisson2", aux)
predC2.novos <- cbind(pred, aux)
## Considering Quasi-Poisson
aux <- exp(confint(
glht(m3QP, linfct = X), calpha = univariate_calpha())$confint)
colnames(aux) <- c("fit", "lwr", "upr")
aux <- data.frame(modelo = "Quasi-Poisson", aux)
predQP.novos <- cbind(pred, aux)
## Representing the confidence intervals
pred.novos <- rbind(predPO.novos, predC1.novos, predC2.novos, predQP.novos)
ord <- order(pred.novos$dose, pred.novos$modelo)
pred.novos <- pred.novos[ord, ]
## Legend
key <- list(columns = 2,
lines = list(col = cols[1:4], lty = rev(1:4), cex = 0.7),
text = list(parse(
text = c("'Poisson'", "'COM-Poisson'",
"'COM-Poisson'[mu]", "'Quasi-Poisson'"))
)
)## Graph
xyplot(novos ~ dose,
data = nitrofen,
xlab = "Nível de adubação potássica",
ylab = "Número de grãos por parcela",
type = c("p", "g"),
alpha = 0.6, key = key,
spread = 0.05,
panel = panel.beeswarm) +
as.layer(
xyplot(fit ~ dose,
auto.key = TRUE,
data = pred.novos,
groups = modelo,
type = "l",
ly = pred.novos$lwr, uy = pred.novos$upr,
cty = "bands", fill = "gray80", alpha = 0.1,
panel = panel.superpose,
panel.groups = panel.cbH,
prepanel = prepanel.cbH,
lty = rev(1:4))
)
## Correlation between estimates
corr.ovos <- purrr::map_df(models.ovos[c("C1", "C2")],
function(x) cov2cor(vcov(x))[1, -1])
corr.ovos## # A tibble: 4 × 2
## C1 C2
## <dbl> <dbl>
## 1 0.99715975 -0.0002671404
## 2 -0.07705505 0.0023319568
## 3 0.15619759 -0.0028685746
## 4 -0.42226644 0.0033258041