##======================================================================
## 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")
}, ...)
}
## 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)
)
## =====================================================================
## 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 = ""))
## =====================================================================
## 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.655
system.time(
m1C2 <- fitcm(form1, data = cottonBolls, model = "CP2", sumto = 50)
)
## user system elapsed
## 3.516 0.000 3.518
m1QP <- 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 39
models.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 ' ' 1
profs.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-05
## =====================================================================
## 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.985
system.time(
m2C2 <- fitcm(form2, data = soyaBeans, model = "CP2", sumto = 700)
)
## user system elapsed
## 14.044 0.004 14.050
m2QP <- 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 72
models.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 ' ' 1
profs.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-05
## =====================================================================
## 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 10
models.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.818
profs.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