|
Applied Statistics to Nematology
|
Meloidogyne Damage Level in Coffee
Walmes Marques Zeviani, Andressa Cristina Zamboni Machado & Santino Aleandro Silva
Descriptive Analysis
#-----------------------------------------------------------------------
# Load packages.
library(lattice)
library(latticeExtra)
library(doBy)
library(multcomp)
library(plyr)
library(wzRfun)
library(nematistics)#-----------------------------------------------------------------------
# A shorter name is better, avoid mispellings.
cp <- coffee_pathogenicity
str(cp)## 'data.frame': 144 obs. of 10 variables:
## $ rootstock: Factor w/ 2 levels "grafted","control": 2 2 2 2 2 2 2 2 2 2 ...
## $ ind : num 0 0 0 0 0 0 0.0625 0.0625 0.0625 0.0625 ...
## $ atw : int 740 680 690 620 810 960 50 220 170 70 ...
## $ frw : num 411 478 415 277 493 ...
## $ vol : int NA NA NA NA NA NA 40 60 60 60 ...
## $ count1 : int NA NA NA NA NA NA 181 1050 1152 469 ...
## $ count2 : int NA NA NA NA NA NA 166 1022 1139 480 ...
## $ fpop : int NA NA NA NA NA NA 13880 124320 137460 56940 ...
## $ nemg : num NA NA NA NA NA ...
## $ rf : num NA NA NA NA NA ...xyplot(atw ~ ind^0.3 | rootstock,
data = cp,
type = c("p", "a"),
grid = TRUE,
xlab = "Initial nematode density (transformed)",
ylab = "Aerial top weight (g)")
xyplot(frw ~ ind^0.3 | rootstock,
data = cp,
type = c("p", "a"),
grid = TRUE,
xlab = "Initial nematode density (transformed)",
ylab = "Fresh root weight (g)")
# Number of nematodes per unit mass (g) of fresh root.
# with(cp, {
# cbind((count1 + count2) * vol/frw, nemg)
# })
xyplot(nemg ~ ind^0.3 | rootstock,
data = cp,
type = c("p", "a"),
grid = TRUE,
scales = list(y = "free"),
xlab = "Initial nematode density (transformed)",
ylab = "Number of nematodes per fresh root unit mass")
#-----------------------------------------------------------------------
# Cell combinations.
xtabs(~ind + rootstock, data = cp)## rootstock
## ind grafted control
## 0 6 6
## 0.0625 6 6
## 0.125 6 6
## 0.25 6 6
## 0.5 6 6
## 1 6 6
## 2 6 6
## 4 6 6
## 8 6 6
## 16 6 6
## 32 6 6
## 64 6 6Aerial Top Weight
#-----------------------------------------------------------------------
# Creating auxiliary variables.
log2(unique(cp$ind))## [1] -Inf -4 -3 -2 -1 0 1 2 3 4 5 6# Replace 0 by 2^{-5} just to plot.
i <- cp$ind == 0
# Indicates if is control (0) or not (1).
cp$ctr <- as.integer(!i)
# Replaces 0 by 2^{-5} to do plots.
cp$ind[i] <- 2^(-5)
# Log base 2 of ind.
cp$lind <- log2(cp$ind)
# Ind as categorical variable.
cp$cind <- factor(cp$lind)#-----------------------------------------------------------------------
# Using ind as a categorical variable.
m0 <- lm(atw ~ rootstock * cind, data = cp)
# Check model assumptions.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# How transform?
MASS::boxcox(m0)
abline(v = 0, col = 2)
# Refit with transformed response.
m1 <- update(m0, log(.) ~ .)
# Check model assumptions again.
par(mfrow = c(2, 2))
plot(m1)
layout(1)
# Anova table.
anova(m1)## Analysis of Variance Table
##
## Response: log(atw)
## Df Sum Sq Mean Sq F value Pr(>F)
## rootstock 1 42.550 42.550 92.6918 < 2.2e-16 ***
## cind 11 40.839 3.713 8.0878 1.958e-10 ***
## rootstock:cind 11 19.785 1.799 3.9182 7.357e-05 ***
## Residuals 120 55.086 0.459
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Hypothesis tests.
lsm <- LSmeans(m1, effect = c("rootstock", "cind"))
lsm## estimate se df t.stat p.value rootstock cind
## 1 6.553177 0.276601 120 23.69181 4.500531e-47 grafted -5
## 2 6.609854 0.276601 120 23.89671 1.915667e-47 control -5
## 3 5.926547 0.276601 120 21.42634 7.891134e-43 grafted -4
## 4 4.550408 0.276601 120 16.45116 1.526840e-32 control -4
## 5 5.788770 0.276601 120 20.92823 7.350186e-42 grafted -3
## 6 4.290312 0.276601 120 15.51083 1.896290e-30 control -3
## 7 5.222639 0.276601 120 18.88149 9.767120e-38 grafted -2
## 8 4.137596 0.276601 120 14.95872 3.373630e-29 control -2
## 9 6.343363 0.276601 120 22.93326 1.108742e-45 grafted -1
## 10 4.285480 0.276601 120 15.49337 2.075999e-30 control -1
## 11 5.349444 0.276601 120 19.33993 1.112381e-38 grafted 0
## 12 4.267314 0.276601 120 15.42769 2.918915e-30 control 0
## 13 5.021688 0.276601 120 18.15499 3.226001e-36 grafted 1
## 14 4.405094 0.276601 120 15.92581 2.229339e-31 control 1
## 15 4.646709 0.276601 120 16.79932 2.630988e-33 grafted 2
## 16 4.434895 0.276601 120 16.03355 1.282974e-31 control 2
## 17 5.225589 0.276601 120 18.89216 9.282943e-38 grafted 3
## 18 4.120550 0.276601 120 14.89709 4.661822e-29 control 3
## 19 5.208805 0.276601 120 18.83148 1.239947e-37 grafted 4
## 20 4.671967 0.276601 120 16.89064 1.662972e-33 control 4
## 21 5.385530 0.276601 120 19.47039 6.023713e-39 grafted 5
## 22 4.588143 0.276601 120 16.58759 7.652261e-33 control 5
## 23 5.917700 0.276601 120 21.39436 9.098643e-43 grafted 6
## 24 3.182271 0.276601 120 11.50492 4.243559e-21 control 6# Get the matrix that computes least squares means.
L <- lsm$K
grid <- lsm$grid
# Split matrix by rootstock levels and update rownames.
L <- by(data = L, INDICES = grid$rootstock, FUN = as.matrix)
L <- lapply(L, FUN = "rownames<-", unique(grid$cind))
# Return the compact letter display for all contrasts.
pmc <- lapply(L, FUN = apmc, model = m1, test = "fdr", focus = "cind")
pmc## $control
## cind fit lwr upr cld
## 1 -5 6.609854 6.062203 7.157504 a
## 2 -4 4.550408 4.002757 5.098059 c
## 3 -3 4.290312 3.742661 4.837962 c
## 4 -2 4.137596 3.589946 4.685247 bc
## 5 -1 4.285480 3.737830 4.833131 c
## 6 0 4.267314 3.719663 4.814965 c
## 7 1 4.405094 3.857443 4.952744 c
## 8 2 4.434895 3.887244 4.982546 c
## 9 3 4.120550 3.572899 4.668201 bc
## 10 4 4.671967 4.124317 5.219618 c
## 11 5 4.588143 4.040492 5.135793 c
## 12 6 3.182271 2.634620 3.729922 b
##
## $grafted
## cind fit lwr upr cld
## 1 -5 6.553177 6.005527 7.100828 b
## 2 -4 5.926547 5.378896 6.474198 bc
## 3 -3 5.788770 5.241120 6.336421 bc
## 4 -2 5.222639 4.674989 5.770290 cd
## 5 -1 6.343363 5.795712 6.891014 ab
## 6 0 5.349444 4.801793 5.897095 acd
## 7 1 5.021688 4.474037 5.569338 cd
## 8 2 4.646709 4.099058 5.194360 d
## 9 3 5.225589 4.677939 5.773240 cd
## 10 4 5.208805 4.661154 5.756456 cd
## 11 5 5.385530 4.837879 5.933180 acd
## 12 6 5.917700 5.370049 6.465350 bc#-----------------------------------------------------------------------
# Stack the list to create a data.frame.
pmc <- ldply(pmc, .id = "rootstock")
str(pmc)## 'data.frame': 24 obs. of 6 variables:
## $ rootstock: Factor w/ 2 levels "control","grafted": 1 1 1 1 1 1 1 1 1 1 ...
## $ cind : Factor w/ 12 levels "0","1","-1","2",..: 11 9 7 5 3 1 2 4 6 8 ...
## $ fit : num 6.61 4.55 4.29 4.14 4.29 ...
## $ lwr : num 6.06 4 3.74 3.59 3.74 ...
## $ upr : num 7.16 5.1 4.84 4.69 4.83 ...
## $ cld : chr "a" "c" "c" "bc" ...# Reorder levels of factors in pmc according to those in cp.
pmc <- equallevels(x = pmc, y = cp)
pmc <- within(pmc, {
txt <- sprintf("%0.2f %s", fit, cld)
xpos <- as.integer(cind) +
0.6 * scale(as.integer(rootstock), scale = FALSE)
})
pmc$cind <- factor(pmc$cind,
labels = c(unique(coffee_pathogenicity$ind)))
key <- list(text = list(levels(pmc$rootstock)),
points = list(pch = c(1, 19)))
segplot(cind ~ lwr + upr,
centers = fit,
data = pmc,
xlab = "Initial nematode density",
ylab = "(log) Aerial top weight",
key = key,
draw = FALSE,
horizontal = FALSE,
groups = rootstock,
gap = 0.1,
pch = key$points$pch[as.integer(pmc$rootstock)],
panel = panel.groups.segplot) +
layer(panel.text(x = pmc$xpos,
y = centers,
labels = pmc$txt,
srt = 90))
Fresh Root Weight
#-----------------------------------------------------------------------
# Using ind as a categorical variable.
m0 <- lm(frw ~ rootstock * cind, data = cp)
# Check model assumptions.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# How transform?
MASS::boxcox(m0)
abline(v = 0, col = 2)
# Refit with transformed response.
m1 <- update(m0, log(.) ~ .)
# Check model assumptions again.
par(mfrow = c(2, 2))
plot(m1)
layout(1)
# Anova table.
anova(m1)## Analysis of Variance Table
##
## Response: log(frw)
## Df Sum Sq Mean Sq F value Pr(>F)
## rootstock 1 58.445 58.445 104.9626 < 2.2e-16 ***
## cind 11 55.842 5.077 9.1170 1.104e-11 ***
## rootstock:cind 11 16.860 1.533 2.7527 0.003254 **
## Residuals 120 66.819 0.557
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Hypothesis tests.
lsm <- LSmeans(m1, effect = c("rootstock", "cind"))
lsm## estimate se df t.stat p.value rootstock cind
## 1 6.298023 0.304637 120 20.673864 2.324662e-41 grafted -5
## 2 6.094644 0.304637 120 20.006253 4.959811e-40 control -5
## 3 5.555450 0.304637 120 18.236298 2.173849e-36 grafted -4
## 4 3.335641 0.304637 120 10.949561 9.086493e-20 control -4
## 5 5.185160 0.304637 120 17.020785 8.663681e-34 grafted -3
## 6 3.445591 0.304637 120 11.310481 1.239636e-20 control -3
## 7 4.154458 0.304637 120 13.637405 3.772967e-26 grafted -2
## 8 3.451299 0.304637 120 11.329219 1.117914e-20 control -2
## 9 5.721312 0.304637 120 18.780756 1.579946e-37 grafted -1
## 10 3.382742 0.304637 120 11.104176 3.869807e-20 control -1
## 11 4.632453 0.304637 120 15.206471 9.231319e-30 grafted 0
## 12 3.522470 0.304637 120 11.562844 3.083989e-21 control 0
## 13 4.347585 0.304637 120 14.271363 1.272160e-27 grafted 1
## 14 3.533482 0.304637 120 11.598994 2.527168e-21 control 1
## 15 4.179688 0.304637 120 13.720225 2.417908e-26 grafted 2
## 16 3.480314 0.304637 120 11.424465 6.611679e-21 control 2
## 17 4.594885 0.304637 120 15.083150 1.758142e-29 grafted 3
## 18 3.249965 0.304637 120 10.668321 4.293914e-19 control 3
## 19 4.478952 0.304637 120 14.702589 1.297213e-28 grafted 4
## 20 3.637056 0.304637 120 11.938984 3.893481e-22 control 4
## 21 4.674479 0.304637 120 15.344425 4.499355e-30 grafted 5
## 22 3.699702 0.304637 120 12.144626 1.259210e-22 control 5
## 23 5.103284 0.304637 120 16.752019 3.338088e-33 grafted 6
## 24 2.802896 0.304637 120 9.200775 1.373562e-15 control 6# Get the matrix that computes least squares means.
L <- lsm$K
grid <- lsm$grid
# Split matrix by rootstock levels and update rownames.
L <- by(data = L, INDICES = grid$rootstock, FUN = as.matrix)
L <- lapply(L, FUN = "rownames<-", unique(grid$cind))
# Return the compact letter display for all contrasts.
pmc <- lapply(L, FUN = apmc, model = m1, test = "fdr", focus = "cind")
pmc## $control
## cind fit lwr upr cld
## 1 -5 6.094644 5.491484 6.697804 a
## 2 -4 3.335641 2.732481 3.938801 b
## 3 -3 3.445591 2.842431 4.048751 b
## 4 -2 3.451299 2.848139 4.054459 b
## 5 -1 3.382742 2.779582 3.985902 b
## 6 0 3.522470 2.919310 4.125630 b
## 7 1 3.533482 2.930322 4.136642 b
## 8 2 3.480314 2.877154 4.083474 b
## 9 3 3.249965 2.646805 3.853125 b
## 10 4 3.637056 3.033896 4.240216 b
## 11 5 3.699702 3.096542 4.302862 b
## 12 6 2.802896 2.199736 3.406056 b
##
## $grafted
## cind fit lwr upr cld
## 1 -5 6.298023 5.694863 6.901183 a
## 2 -4 5.555450 4.952290 6.158610 ac
## 3 -3 5.185160 4.582000 5.788320 bcd
## 4 -2 4.154458 3.551298 4.757618 d
## 5 -1 5.721312 5.118152 6.324472 ab
## 6 0 4.632453 4.029293 5.235613 cd
## 7 1 4.347585 3.744425 4.950745 d
## 8 2 4.179688 3.576528 4.782847 d
## 9 3 4.594885 3.991725 5.198045 cd
## 10 4 4.478952 3.875792 5.082112 d
## 11 5 4.674479 4.071319 5.277639 bcd
## 12 6 5.103284 4.500124 5.706444 bcd#-----------------------------------------------------------------------
# Stack the list to create a data.frame.
pmc <- ldply(pmc, .id = "rootstock")
str(pmc)## 'data.frame': 24 obs. of 6 variables:
## $ rootstock: Factor w/ 2 levels "control","grafted": 1 1 1 1 1 1 1 1 1 1 ...
## $ cind : Factor w/ 12 levels "0","1","-1","2",..: 11 9 7 5 3 1 2 4 6 8 ...
## $ fit : num 6.09 3.34 3.45 3.45 3.38 ...
## $ lwr : num 5.49 2.73 2.84 2.85 2.78 ...
## $ upr : num 6.7 3.94 4.05 4.05 3.99 ...
## $ cld : chr "a" "b" "b" "b" ...# Reorder levels of factors in pmc according to those in cp.
pmc <- equallevels(x = pmc, y = cp)
pmc <- within(pmc, {
txt <- sprintf("%0.2f %s", fit, cld)
xpos <- as.integer(cind) +
0.6 * scale(as.integer(rootstock), scale = FALSE)
})
pmc$cind <- factor(pmc$cind,
labels = c(unique(coffee_pathogenicity$ind)))
key <- list(text = list(levels(pmc$rootstock)),
points = list(pch = c(1, 19)))
segplot(cind ~ lwr + upr,
centers = fit,
data = pmc,
xlab = "Initial nematode density",
ylab = "(log) Fresh root weight",
key = key,
draw = FALSE,
horizontal = FALSE,
groups = rootstock,
gap = 0.1,
pch = key$points$pch[as.integer(pmc$rootstock)],
panel = panel.groups.segplot) +
layer(panel.text(x = pmc$xpos,
y = centers,
labels = pmc$txt,
srt = 90))
Nematode Density
Analysis with Gaussian linear model
#-----------------------------------------------------------------------
# Dropout lines with NA.
cp <- droplevels(na.omit(cp))
m0 <- lm(nemg ~ rootstock * cind, data = cp)
# Check model assumptions.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# How transform?
MASS::boxcox(m0)
abline(v = 0, col = 2)
# Refit with transformed response.
m1 <- update(m0, log(.) ~ .)
# Check model assumptions again.
par(mfrow = c(2, 2))
plot(m1)
layout(1)
# Anova table.
anova(m1)## Analysis of Variance Table
##
## Response: log(nemg)
## Df Sum Sq Mean Sq F value Pr(>F)
## rootstock 1 267.652 267.652 344.969 < 2.2e-16 ***
## cind 10 72.009 7.201 9.281 5.301e-11 ***
## rootstock:cind 10 122.583 12.258 15.799 < 2.2e-16 ***
## Residuals 110 85.346 0.776
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Hypothesis tests.
lsm <- LSmeans(m1, effect = c("rootstock", "cind"))
lsm## estimate se df t.stat p.value rootstock cind
## 1 1.593286 0.3595998 110 4.430720 2.229201e-05 grafted -4
## 2 7.263462 0.3595998 110 20.198741 8.317915e-39 control -4
## 3 2.365406 0.3595998 110 6.577885 1.668063e-09 grafted -3
## 4 4.854279 0.3595998 110 13.499115 4.390882e-25 control -3
## 5 4.345322 0.3595998 110 12.083773 6.623302e-22 grafted -2
## 6 3.684867 0.3595998 110 10.247135 1.056739e-17 control -2
## 7 2.624409 0.3595998 110 7.298140 4.826215e-11 grafted -1
## 8 3.093329 0.3595998 110 8.602144 6.088220e-14 control -1
## 9 1.961710 0.3595998 110 5.455259 3.043494e-07 grafted 0
## 10 7.580697 0.3595998 110 21.080930 1.965192e-40 control 0
## 11 2.272724 0.3595998 110 6.320149 5.725926e-09 grafted 1
## 12 4.486912 0.3595998 110 12.477517 8.505880e-23 control 1
## 13 2.240829 0.3595998 110 6.231454 8.711582e-09 grafted 2
## 14 7.227265 0.3595998 110 20.098082 1.283071e-38 control 2
## 15 3.138345 0.3595998 110 8.727328 3.169109e-14 grafted 3
## 16 6.641504 0.3595998 110 18.469157 1.701315e-35 control 3
## 17 2.529988 0.3595998 110 7.035566 1.782314e-10 grafted 4
## 18 4.967106 0.3595998 110 13.812874 8.883318e-26 control 4
## 19 3.696855 0.3595998 110 10.280472 8.858626e-18 grafted 5
## 20 6.684461 0.3595998 110 18.588614 9.928039e-36 control 5
## 21 2.397603 0.3595998 110 6.667420 1.081792e-09 grafted 6
## 22 4.009733 0.3595998 110 11.150544 8.914645e-20 control 6# Get the matrix that computes least squares means.
L <- lsm$K
grid <- lsm$grid
# Split matrix by rootstock levels and update rownames.
L <- by(data = L, INDICES = grid$rootstock, FUN = as.matrix)
L <- lapply(L, FUN = "rownames<-", unique(grid$cind))
# Return the compact letter display for all contrasts.
pmc <- lapply(L, FUN = apmc, model = m1, test = "fdr", focus = "cind")
pmc## $control
## cind fit lwr upr cld
## 1 -4 7.263462 6.550820 7.976105 a
## 2 -3 4.854279 4.141636 5.566921 c
## 3 -2 3.684867 2.972225 4.397509 bd
## 4 -1 3.093329 2.380686 3.805971 d
## 5 0 7.580697 6.868055 8.293339 a
## 6 1 4.486912 3.774270 5.199554 bc
## 7 2 7.227265 6.514623 7.939907 a
## 8 3 6.641504 5.928862 7.354146 a
## 9 4 4.967106 4.254464 5.679748 c
## 10 5 6.684461 5.971819 7.397103 a
## 11 6 4.009733 3.297091 4.722375 cd
##
## $grafted
## cind fit lwr upr cld
## 1 -4 1.593286 0.8806434 2.305928 c
## 2 -3 2.365406 1.6527635 3.078048 cd
## 3 -2 4.345322 3.6326795 5.057964 a
## 4 -1 2.624409 1.9117669 3.337051 bcd
## 5 0 1.961710 1.2490674 2.674352 cd
## 6 1 2.272724 1.5600818 2.985366 cd
## 7 2 2.240829 1.5281870 2.953472 cd
## 8 3 3.138345 2.4257027 3.850987 ad
## 9 4 2.529988 1.8173456 3.242630 bcd
## 10 5 3.696855 2.9842128 4.409497 ab
## 11 6 2.397603 1.6849602 3.110245 cd#-----------------------------------------------------------------------
# Stack the list to create a data.frame.
pmc <- ldply(pmc, .id = "rootstock")
str(pmc)## 'data.frame': 22 obs. of 6 variables:
## $ rootstock: Factor w/ 2 levels "control","grafted": 1 1 1 1 1 1 1 1 1 1 ...
## $ cind : Factor w/ 11 levels "0","1","-1","2",..: 9 7 5 3 1 2 4 6 8 10 ...
## $ fit : num 7.26 4.85 3.68 3.09 7.58 ...
## $ lwr : num 6.55 4.14 2.97 2.38 6.87 ...
## $ upr : num 7.98 5.57 4.4 3.81 8.29 ...
## $ cld : chr "a" "c" "bd" "d" ...# Reorder levels of factors in pmc according to those in cp.
pmc <- droplevels(equallevels(x = pmc, y = cp))
pmc <- within(pmc, {
txt <- sprintf("%0.2f %s", fit, cld)
xpos <- as.integer(cind) +
0.6 * scale(as.integer(rootstock), scale = FALSE)
})
pmc$cind <- factor(pmc$cind,
labels = c(unique(coffee_pathogenicity$ind)[-1]))
key <- list(text = list(levels(pmc$rootstock)),
points = list(pch = c(1, 19)))
segplot(cind ~ lwr + upr,
centers = fit,
data = pmc,
xlab = "Initial nematode density",
ylab = "(log) Nematode by fresh root unit mass",
key = key,
draw = FALSE,
horizontal = FALSE,
groups = rootstock,
gap = 0.1,
pch = key$points$pch[as.integer(pmc$rootstock)],
panel = panel.groups.segplot) +
layer(panel.text(x = pmc$xpos,
y = centers,
labels = pmc$txt,
srt = 90))
Analysis with Quasi-likelihood log-linear model
#-----------------------------------------------------------------------
cp$count <- with(cp, count1 + count2)
cp$off <- with(cp, frw/vol)
str(cp)## 'data.frame': 132 obs. of 15 variables:
## $ rootstock: Factor w/ 2 levels "grafted","control": 2 2 2 2 2 2 2 2 2 2 ...
## $ ind : num 0.0625 0.0625 0.0625 0.0625 0.0625 0.0625 0.125 0.125 0.125 0.125 ...
## $ atw : int 50 220 170 70 110 50 50 70 60 50 ...
## $ frw : num 15 40.6 72.7 19.8 32.8 ...
## $ vol : int 40 60 60 60 60 60 60 60 60 60 ...
## $ count1 : int 181 1050 1152 469 127 186 45 36 22 41 ...
## $ count2 : int 166 1022 1139 480 98 194 51 29 31 36 ...
## $ fpop : int 13880 124320 137460 56940 13500 22800 5760 3900 3180 4620 ...
## $ nemg : num 925 3061 1891 2870 411 ...
## $ rf : num 11.1 99.5 110 45.6 10.8 ...
## $ ctr : int 1 1 1 1 1 1 1 1 1 1 ...
## $ lind : num -4 -4 -4 -4 -4 -4 -3 -3 -3 -3 ...
## $ cind : Factor w/ 11 levels "-4","-3","-2",..: 1 1 1 1 1 1 2 2 2 2 ...
## $ count : int 347 2072 2291 949 225 380 96 65 53 77 ...
## $ off : num 0.375 0.677 1.212 0.331 0.547 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:12] 1 2 3 4 5 6 73 74 75 76 ...
## .. ..- attr(*, "names")= chr [1:12] "1" "2" "3" "4" ...m0 <- glm(count ~ offset(log(off)) + rootstock * cind,
data = cp,
family = quasi(variance = "mu^2", link = "log"))
m1 <- glm(count ~ log(off) + rootstock * cind,
data = cp,
family = quasi(variance = "mu^2", link = "log"))
anova(m0, m1, test = "F")## Analysis of Deviance Table
##
## Model 1: count ~ offset(log(off)) + rootstock * cind
## Model 2: count ~ log(off) + rootstock * cind
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 110 70.181
## 2 109 25.233 1 44.948 222.93 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Check model assumptions.
par(mfrow = c(2, 2))
plot(m1)
layout(1)
summary(m1)##
## Call:
## glm(formula = count ~ log(off) + rootstock * cind, family = quasi(variance = "mu^2",
## link = "log"), data = cp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.21573 -0.35464 -0.02664 0.25037 0.85321
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.61252 0.19088 13.687 < 2e-16 ***
## log(off) 0.06664 0.05599 1.190 0.236524
## rootstockcontrol 4.36413 0.27505 15.867 < 2e-16 ***
## cind-3 0.43234 0.26007 1.662 0.099313 .
## cind-2 1.52980 0.27085 5.648 1.31e-07 ***
## cind-1 1.03054 0.25931 3.974 0.000127 ***
## cind0 -0.63532 0.26503 -2.397 0.018224 *
## cind1 -0.12196 0.26332 -0.463 0.644160
## cind2 -0.74848 0.27045 -2.768 0.006637 **
## cind3 0.58395 0.26476 2.206 0.029515 *
## cind4 -0.05998 0.26616 -0.225 0.822137
## cind5 1.15964 0.26390 4.394 2.59e-05 ***
## cind6 0.39097 0.26048 1.501 0.136259
## rootstockcontrol:cind-3 -3.13168 0.36736 -8.525 9.63e-14 ***
## rootstockcontrol:cind-2 -5.14633 0.37896 -13.580 < 2e-16 ***
## rootstockcontrol:cind-1 -5.55258 0.36695 -15.132 < 2e-16 ***
## rootstockcontrol:cind0 0.91564 0.37362 2.451 0.015848 *
## rootstockcontrol:cind1 -2.67382 0.37211 -7.186 8.76e-11 ***
## rootstockcontrol:cind2 0.63010 0.37789 1.667 0.098298 .
## rootstockcontrol:cind3 -1.53715 0.37006 -4.154 6.52e-05 ***
## rootstockcontrol:cind4 -2.13081 0.37603 -5.667 1.20e-07 ***
## rootstockcontrol:cind5 -1.50386 0.37411 -4.020 0.000108 ***
## rootstockcontrol:cind6 -4.37330 0.36663 -11.928 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasi family taken to be 0.2016245)
##
## Null deviance: 406.727 on 131 degrees of freedom
## Residual deviance: 25.233 on 109 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 7# Testa se o \beta do offset é 1.
L <- rep(0, length(coef(m1)))
L[2] <- 1
summary(glht(m1, linfct = rbind(L), rhs = 1))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: glm(formula = count ~ log(off) + rootstock * cind, family = quasi(variance = "mu^2",
## link = "log"), data = cp)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## L == 1 0.06664 0.05599 -16.67 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# m1 <- m0
# Anova table.
anova(m1, test = "F")## Analysis of Deviance Table
##
## Model: quasi, link: log
##
## Response: count
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 131 406.73
## log(off) 1 16.617 130 390.11 82.418 5.394e-15
## rootstock 1 197.567 129 192.54 979.875 < 2.2e-16
## cind 10 54.230 119 138.31 26.896 < 2.2e-16
## rootstock:cind 10 113.080 109 25.23 56.085 < 2.2e-16
##
## NULL
## log(off) ***
## rootstock ***
## cind ***
## rootstock:cind ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Hypothesis tests.
# Get the matrix that computes least squares means.
grid <- unique(subset(cp, select = c("rootstock", "cind")))
L <- model.matrix(formula(m1)[-2],
data = cbind(off = 1, grid))
# Split matrix by rootstock levels and update rownames.
L <- by(data = L, INDICES = grid$rootstock, FUN = as.matrix)
L <- lapply(L, FUN = "rownames<-", unique(grid$cind))
# Return the compact letter display for all contrasts.
pmc <- lapply(L, FUN = apmc, model = m1, test = "fdr", focus = "cind")
pmc## $grafted
## cind fit lwr upr cld
## 1 -4 2.612522 2.238407 2.986637 bg
## 2 -3 3.044858 2.679976 3.409741 bf
## 3 -2 4.142319 3.779642 4.504997 e
## 4 -1 3.643058 3.265567 4.020550 de
## 5 0 1.977200 1.617892 2.336508 ac
## 6 1 2.490559 2.131003 2.850114 cg
## 7 2 1.864040 1.501730 2.226350 a
## 8 3 3.196474 2.837183 3.555765 df
## 9 4 2.552546 2.192990 2.912102 bg
## 10 5 3.772161 3.412791 4.131530 e
## 11 6 3.003488 2.640065 3.366911 fg
##
## $control
## cind fit lwr upr cld
## 1 -4 6.976654 6.609448 7.343860 ab
## 2 -3 4.277309 3.911035 4.643582 de
## 3 -2 3.360122 2.998427 3.721817 g
## 4 -1 2.454606 2.083178 2.826034 f
## 5 0 7.256973 6.895150 7.618796 a
## 6 1 4.180868 3.819185 4.542550 e
## 7 2 6.858272 6.495873 7.220670 ab
## 8 3 6.023458 5.654785 6.392131 c
## 9 4 4.785864 4.425310 5.146418 d
## 10 5 6.632433 6.272158 6.992708 b
## 11 6 2.994316 2.614973 3.373658 g#-----------------------------------------------------------------------
# Stack the list to create a data.frame.
pmc <- ldply(pmc, .id = "rootstock")
str(pmc)## 'data.frame': 22 obs. of 6 variables:
## $ rootstock: Factor w/ 2 levels "grafted","control": 1 1 1 1 1 1 1 1 1 1 ...
## $ cind : Factor w/ 11 levels "0","1","-1","2",..: 9 7 5 3 1 2 4 6 8 10 ...
## $ fit : num 2.61 3.04 4.14 3.64 1.98 ...
## $ lwr : num 2.24 2.68 3.78 3.27 1.62 ...
## $ upr : num 2.99 3.41 4.5 4.02 2.34 ...
## $ cld : chr "bg" "bf" "e" "de" ...# Reorder levels of factors in pmc according to those in cp.
pmc <- droplevels(equallevels(x = pmc, y = cp))
pmc <- within(pmc, {
txt <- sprintf("%0.2f %s", fit, cld)
xpos <- as.integer(cind) +
0.6 * scale(as.integer(rootstock), scale = FALSE)
})
key <- list(text = list(levels(pmc$rootstock)),
points = list(pch = c(1, 19)))
segplot(cind ~ lwr + upr,
centers = fit,
data = pmc,
xlab = "Initial nematode density",
ylab = "(log) Nematode by fresh root unit mass",
key = key,
draw = FALSE,
horizontal = FALSE,
groups = rootstock,
gap = 0.1,
pch = key$points$pch[as.integer(pmc$rootstock)],
panel = panel.groups.segplot) +
layer(panel.text(x = pmc$xpos,
y = centers,
labels = pmc$txt,
srt = 90))
Session Information
## Updated in 2016-08-15.
##
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.1 LTS
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods
## [7] base
##
## other attached packages:
## [1] plyr_1.8.4 multcomp_1.4-6 TH.data_1.0-7
## [4] MASS_7.3-45 survival_2.39-4 mvtnorm_1.0-5
## [7] doBy_4.5-15 latticeExtra_0.6-28 RColorBrewer_1.1-2
## [10] lattice_0.20-33 knitr_1.13 nematistics_0.0-1
## [13] wzRfun_0.70 gsubfn_0.6-6 proto_0.3-10
## [16] devtools_1.11.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.3 formatR_1.4 git2r_0.15.0
## [4] tools_3.3.1 digest_0.6.9 memoise_1.0.0
## [7] evaluate_0.9 Matrix_1.2-6 curl_0.9.7
## [10] yaml_2.1.13 rpanel_1.1-3 withr_1.0.1
## [13] stringr_1.0.0 httr_1.2.0 roxygen2_5.0.1
## [16] grid_3.3.1 R6_2.1.2 tcltk_3.3.1
## [19] rmarkdown_0.9.6 magrittr_1.5 codetools_0.2-14
## [22] htmltools_0.3.5 splines_3.3.1 sandwich_2.3-4
## [25] stringi_1.1.1 zoo_1.7-13
Applied Statistics to Nematology |
github.com/walmes/nematistics |