|
Multivariate Covariance Generalized Linear Models for the Analysis of Experimental Data
|
Water and potassium on the production of soybean
Wagner H. Bonat
Walmes M. Zeviani
wbonat@ufpr.brwalmes@ufpr.br
LEG/DEST/UFPR
LEG/DEST/UFPR
wbonat@ufpr.brwalmes@ufpr.br62a RBras & 17o SEAGRO July 24–28, 2017 UFLA, Lavras/MG
1 Data description and objectives
The tropical soils, usually poor in potassium (K), demand potassium fertilization when cultivated with soybean (Glycine max L.) to obtain satisfactory yields. The aim of this study was to evaluate the effects of K doses (pts) and soil humidity (wtr) levels on soybean agronomic characteristics.
The experiment was carried out in a greenhouse, in pots with two plants, containing 5 dm3 of soil. The experimental design was completely randomized block with treatments in a 5 \(\times\) 3 factorial arrangement. The K doses were 0; 30; 60, 120 and 180 mg dm-3 and the soil humidity ranged from 35 to 40; 47.5 to 52.5; and 60 to 65% of the total porosity. The responses measured were: grain yield (yield), weight of a hundred grains (w100), total number of grains per pot (tg), K level in the grain (Kconc), number of viable pods (nvp, nip is the inviable pods).
For more information, see Serafim et al. (2012).
#-----------------------------------------------------------------------
# Packages.
rm(list = ls())
library(lattice)
library(latticeExtra)
library(car)
library(candisc)
library(corrplot)
library(doBy)
library(multcomp)
library(mcglm)
library(Matrix)2 MLM for the soybean dataset
#-----------------------------------------------------------------------
# Multivariate analysis of the soybean dataset.
# Visit to see the online documentation:
# https://github.com/walmes/wzRfun/blob/master/R/wzRfun.R#L10.
data(soybeanwp, package = "wzRfun")
str(soybeanwp)## 'data.frame': 75 obs. of 9 variables:
## $ potassium: int 0 30 60 120 180 0 30 60 120 180 ...
## $ water : num 37.5 37.5 37.5 37.5 37.5 50 50 50 50 50 ...
## $ block : Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ yield : num 14.6 21.5 24.6 21.9 28.1 ...
## $ w100 : num 10.7 13.5 15.8 12.8 14.8 ...
## $ Kconc : num 15.1 17.1 19.1 18.1 19.1 ...
## $ tg : int 136 159 156 171 190 140 193 200 208 237 ...
## $ nip : int 22 2 0 2 0 20 6 6 7 10 ...
## $ nvp : int 56 62 66 68 82 63 86 94 86 97 ...head(soybeanwp)## potassium water block yield w100 Kconc tg nip nvp
## 1 0 37.5 I 14.55 10.70 15.13 136 22 56
## 2 30 37.5 I 21.51 13.53 17.12 159 2 62
## 3 60 37.5 I 24.62 15.78 19.11 156 0 66
## 4 120 37.5 I 21.88 12.80 18.12 171 2 68
## 5 180 37.5 I 28.11 14.79 19.11 190 0 82
## 6 0 50.0 I 17.16 12.26 12.14 140 20 63# A copy with shorter names.
swp <- soybeanwp
# dput(names(swp))
names(swp) <- c("pts", "wtr", "blk",
"yield", "w100", "Kconc", "tg", "nip", "nvp")
# Convert the controled numeric factors to categorical factors.
swp <- transform(swp,
Pts = factor(pts),
Wtr = factor(wtr))
# Creates the proportion of viable pods.
swp$pvp <- with(swp, nvp/(nvp + nip))
str(swp)## 'data.frame': 75 obs. of 12 variables:
## $ pts : int 0 30 60 120 180 0 30 60 120 180 ...
## $ wtr : num 37.5 37.5 37.5 37.5 37.5 50 50 50 50 50 ...
## $ blk : Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ yield: num 14.6 21.5 24.6 21.9 28.1 ...
## $ w100 : num 10.7 13.5 15.8 12.8 14.8 ...
## $ Kconc: num 15.1 17.1 19.1 18.1 19.1 ...
## $ tg : int 136 159 156 171 190 140 193 200 208 237 ...
## $ nip : int 22 2 0 2 0 20 6 6 7 10 ...
## $ nvp : int 56 62 66 68 82 63 86 94 86 97 ...
## $ Pts : Factor w/ 5 levels "0","30","60",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ Wtr : Factor w/ 3 levels "37.5","50","62.5": 1 1 1 1 1 2 2 2 2 2 ...
## $ pvp : num 0.718 0.969 1 0.971 1 ...#-----------------------------------------------------------------------
# Data visulization.
# The 74 observation is an outlier to yield and tg.
combineLimits(
useOuterStrips(
xyplot(yield + tg + w100 + Kconc + pvp ~ pts | Wtr,
outer = TRUE,
as.table = TRUE,
# data = swp[-74, ],
data = swp,
pch = 19,
xlab = "Potassium",
ylab = "Values",
scales = list(y = "free"),
groups = blk))) +
layer(panel.xyplot(x = x, y = y, type = "a", col = 1))
# Data summaries.
summary(swp)## pts wtr blk yield w100
## Min. : 0 Min. :37.5 I :15 Min. :10.30 Min. :10.10
## 1st Qu.: 30 1st Qu.:37.5 II :15 1st Qu.:19.71 1st Qu.:12.25
## Median : 60 Median :50.0 III:15 Median :24.44 Median :13.34
## Mean : 78 Mean :50.0 IV :15 Mean :24.77 Mean :13.69
## 3rd Qu.:120 3rd Qu.:62.5 V :15 3rd Qu.:29.78 3rd Qu.:15.12
## Max. :180 Max. :62.5 Max. :41.04 Max. :18.23
## Kconc tg nip nvp
## Min. :10.65 Min. : 92.0 Min. : 0.000 Min. : 36.00
## 1st Qu.:15.13 1st Qu.:147.0 1st Qu.: 1.000 1st Qu.: 61.50
## Median :18.61 Median :170.0 Median : 3.000 Median : 71.00
## Mean :17.79 Mean :180.3 Mean : 6.093 Mean : 75.23
## 3rd Qu.:20.36 3rd Qu.:215.5 3rd Qu.: 8.000 3rd Qu.: 90.50
## Max. :23.10 Max. :271.0 Max. :27.000 Max. :110.00
## Pts Wtr pvp
## 0 :15 37.5:25 Min. :0.6538
## 30 :15 50 :25 1st Qu.:0.9106
## 60 :15 62.5:25 Median :0.9636
## 120:15 Mean :0.9234
## 180:15 3rd Qu.:0.9848
## Max. :1.00002.1 Analysing the subset for potassium equals 0
#-----------------------------------------------------------------------
# Condition 1: no sufficient evidence in each response.
# Subsetting the dataset.
da <- subset(swp[-74, ], pts == 0)
xtabs(~blk + Wtr, data = da)## Wtr
## blk 37.5 50 62.5
## I 1 1 1
## II 1 1 1
## III 1 1 1
## IV 1 1 1
## V 1 1 1xyplot(yield + tg + w100 + Kconc + pvp ~ wtr,
outer = TRUE,
scales = "free",
type = c("p", "a"),
data = da)
splom(~da[c(4:7, 12)],
data = da,
groups = Wtr,
as.matrix = TRUE)
# `Wtr` is a 3 level categorical factor.
m0 <- lm(cbind(yield, tg, w100, Kconc, pvp) ~ blk + Wtr,
data = da)
# Separated anovas.
summary.aov(m0)## Response yield :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 4.7999 1.2000 0.3937 0.80798
## Wtr 2 19.2544 9.6272 3.1585 0.09749 .
## Residuals 8 24.3842 3.0480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response tg :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 456.4 114.10 0.7591 0.57990
## Wtr 2 1462.9 731.47 4.8667 0.04142 *
## Residuals 8 1202.4 150.30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response w100 :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 2.6981 0.67452 1.3088 0.3447
## Wtr 2 0.2942 0.14709 0.2854 0.7591
## Residuals 8 4.1230 0.51538
##
## Response Kconc :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 2.5354 0.6339 0.2837 0.88053
## Wtr 2 14.5229 7.2615 3.2506 0.09263 .
## Residuals 8 17.8711 2.2339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response pvp :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 0.006158 0.0015394 0.2018 0.9304
## Wtr 2 0.017134 0.0085670 1.1229 0.3717
## Residuals 8 0.061033 0.0076291# Multivariate anova.
Anova(m0)##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 1.23325 0.62403 20 28 0.861
## Wtr 2 0.99592 0.99187 10 10 0.505Anova(m0, test = "Roy")##
## Type II MANOVA Tests: Roy test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 5.2653 7.3714 5 7 0.01034 *
## Wtr 2 3.8746 3.8746 5 5 0.08173 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Residuals.
# Raw residuals.
r <- residuals(m0)
cor(r)## yield tg w100 Kconc pvp
## yield 1.0000000 0.84447670 0.60857389 -0.8161456 -0.2705549
## tg 0.8444767 1.00000000 0.09008227 -0.7158224 -0.2494919
## w100 0.6085739 0.09008227 1.00000000 -0.4753493 -0.1545910
## Kconc -0.8161456 -0.71582244 -0.47534932 1.0000000 0.2009581
## pvp -0.2705549 -0.24949193 -0.15459096 0.2009581 1.0000000# Change the 3rd color of the palette used for the ellipses.
oldpal <- palette()
palette(c("black", "blue", "red"))
scatterplotMatrix(r,
gap = 0,
smooth = FALSE,
reg.line = FALSE,
ellipse = TRUE,
diagonal = "qqplot")
palette(oldpal)
corrplot(cor(r),
type = "upper",
tl.pos = "d",
outline = TRUE,
method = "ellipse")
#-----------------------------------------------------------------------
# `wtr` is a 3 level numerical factor.
m1 <- update(m0, . ~ blk + wtr + I(wtr^2))
# Separated anovas.
summary.aov(m1)## Response yield :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 4.7999 1.2000 0.3937 0.80798
## wtr 1 16.5123 16.5123 5.4174 0.04835 *
## I(wtr^2) 1 2.7422 2.7422 0.8997 0.37064
## Residuals 8 24.3842 3.0480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response tg :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 456.40 114.10 0.7591 0.57990
## wtr 1 1392.40 1392.40 9.2641 0.01597 *
## I(wtr^2) 1 70.53 70.53 0.4693 0.51267
## Residuals 8 1202.40 150.30
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response w100 :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 2.6981 0.67452 1.3088 0.3447
## wtr 1 0.0548 0.05476 0.1063 0.7528
## I(wtr^2) 1 0.2394 0.23941 0.4645 0.5148
## Residuals 8 4.1230 0.51538
##
## Response Kconc :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 2.5354 0.6339 0.2837 0.88053
## wtr 1 13.1102 13.1102 5.8688 0.04168 *
## I(wtr^2) 1 1.4127 1.4127 0.6324 0.44943
## Residuals 8 17.8711 2.2339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response pvp :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 0.006158 0.0015394 0.2018 0.9304
## wtr 1 0.005078 0.0050776 0.6656 0.4382
## I(wtr^2) 1 0.012056 0.0120563 1.5803 0.2442
## Residuals 8 0.061033 0.0076291# Multivariate anova.
anova(m1) # Sequential (more appropriate for the current case).## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 1.00000 231498 5 4 5.225e-11 ***
## blk 4 1.23325 1 20 28 0.8610
## wtr 1 0.76726 3 5 4 0.1843
## I(wtr^2) 1 0.45346 1 5 4 0.6723
## Residuals 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(m1) # Marginal.##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 1.23325 0.62403 20 28 0.8610
## wtr 1 0.50811 0.82638 5 4 0.5897
## I(wtr^2) 1 0.45346 0.66376 5 4 0.6723anova(m1, test = "Roy") # Sequential (more appropriate ...).## Analysis of Variance Table
##
## Df Roy approx F num Df den Df Pr(>F)
## (Intercept) 1 289366 231493 5 4 5.225e-11 ***
## blk 4 5 7 5 7 0.01034 *
## wtr 1 3 3 5 4 0.18432
## I(wtr^2) 1 1 1 5 4 0.67233
## Residuals 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(m1, test = "Roy") # Marginal.##
## Type II MANOVA Tests: Roy test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 5.2653 7.3714 5 7 0.01034 *
## wtr 1 1.0330 0.8264 5 4 0.58966
## I(wtr^2) 1 0.8297 0.6638 5 4 0.67233
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# NOTE: Univariate and Multivariate in disagreement! The way of
# representing the factors (numerical or categorical) is important. BUT
# in this case, the difference is due the type of hypothesis display
# (type I, II or III).
#-----------------------------------------------------------------------
# Performing univariate analysis.
tb <- lapply(c("yield", "tg", "w100", "Kconc", "pvp"),
FUN = function(y) {
m0 <- lm(as.formula(sprintf("%s ~ blk + Wtr", y)),
data = da)
m <- LSmeans(m0,
effect = "Wtr")
g <- glht(model = m0,
linfct = mcp(Wtr = "Tukey"))
g <- cld(g, decreasing = TRUE)$mcletters$Letters
r <- data.frame(
Wtr = m$grid$Wtr,
means = sprintf("%0.2f %s",
m$coef[, "estimate"],
g[as.character(m$grid$Wtr)]))
names(r)[2] <- y
return(r)
})
# Table of means (recursive merging to put the results together).
tbm <- Reduce(merge, x = tb)
tbm## Wtr yield tg w100 Kconc pvp
## 1 37.5 13.52 a 115.80 b 11.69 a 15.03 a 0.81 a
## 2 50 15.71 a 132.20 ab 11.89 a 13.23 a 0.73 a
## 3 62.5 16.09 a 139.40 a 11.54 a 12.74 a 0.76 a#-----------------------------------------------------------------------
# Canonical discriminant analysis.
cd_W <- candisc(m0, term = "Wtr")
cd_W##
## Canonical Discriminant Analysis for Wtr:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.79485 3.87460 3.6229 93.9009 93.901
## 2 0.20107 0.25167 3.6229 6.0991 100.000
##
## Test of H0: The canonical correlations in the
## current row and all that follow are zero
##
## LR test stat approx F numDF denDF Pr(> F)
## 1 0.16390 2.35215 10 16 0.06149 .
## 2 0.79893 0.56625 4 9 0.69366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(cd_W)##
## Canonical Discriminant Analysis for Wtr:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.79485 3.87460 3.6229 93.9009 93.901
## 2 0.20107 0.25167 3.6229 6.0991 100.000
##
## Class means:
##
## Can1 Can2
## 37.5 2.02006 -0.058267
## 50 -0.81203 0.474989
## 62.5 -1.20802 -0.416722
##
## std coefficients:
## Can1 Can2
## yield -33.91027 11.325694
## tg 27.43927 -9.714454
## w100 18.86256 -5.638358
## Kconc 1.27661 -0.061395
## pvp 0.55525 -0.801308# The weights to compute the scores.
cd_W$coeffs.raw## Can1 Can2
## yield -19.4232570 6.48717611
## tg 2.2381701 -0.79238989
## w100 26.2746880 -7.85397650
## Kconc 0.8541409 -0.04107708
## pvp 6.3569486 -9.17405677# The scores.
head(cd_W$scores)## blk Wtr Can1 Can2
## 1 I 37.5 0.6692268 -0.76325088
## 6 I 50 -2.3769713 0.49239795
## 11 I 62.5 -1.8374999 -0.21913727
## 16 II 37.5 3.2878932 1.03290374
## 21 II 50 2.2521080 -0.05178607
## 26 II 62.5 -0.1789190 -1.58288545# Calculating the scores (with the centered y-variables).
Y <- m0$model[[1]]
z <- scale(Y, scale = FALSE) %*% cd_W$coeffs.raw
head(z[, 1:2])## Can1 Can2
## [1,] 0.6692268 -0.76325088
## [2,] -2.3769713 0.49239795
## [3,] -1.8374999 -0.21913727
## [4,] 3.2878932 1.03290374
## [5,] 2.2521080 -0.05178607
## [6,] -0.1789190 -1.58288545head(cd_W$scores[, c("Can1", "Can2")])## Can1 Can2
## 1 0.6692268 -0.76325088
## 6 -2.3769713 0.49239795
## 11 -1.8374999 -0.21913727
## 16 3.2878932 1.03290374
## 21 2.2521080 -0.05178607
## 26 -0.1789190 -1.58288545# The loadings are the correlations between scores and responses.
cd_W$structure## Can1 Can2
## yield -0.4836916 -0.07599147
## tg -0.7093883 -0.33842087
## w100 0.3986073 0.45934878
## Kconc 0.6111820 0.17628239
## pvp 0.3865570 -0.50617860cor(cbind(Y, z))[colnames(Y), colnames(z)]## Can1 Can2
## yield -0.4836916 -0.07599147
## tg -0.7093883 -0.33842087
## w100 0.3986073 0.45934878
## Kconc 0.6111820 0.17628239
## pvp 0.3865570 -0.50617860# The biplot.
plot(cd_W)## Vector scale factor set to 4.734
2.2 Analysing the subset for potassium equals 120
#-----------------------------------------------------------------------
# Condition 2: results are similar for sets of responses.
# Subsetting the dataset.
da <- subset(swp[-74, ], pts == 120)
xtabs(~blk + Wtr, data = da)## Wtr
## blk 37.5 50 62.5
## I 1 1 1
## II 1 1 1
## III 1 1 1
## IV 1 1 1
## V 1 1 0xyplot(yield + tg + w100 + Kconc + pvp ~ wtr,
outer = TRUE,
scales = "free",
type = c("p", "a"),
data = da)
splom(~da[c(4:7, 12)],
data = da,
groups = Wtr,
as.matrix = TRUE)
# `Wtr` is a 3 level categorical factor.
m0 <- lm(cbind(yield, tg, w100, Kconc, pvp) ~ blk + Wtr,
data = da)
# Separated anovas.
summary.aov(m0)## Response yield :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 7.247 1.812 0.2282 0.914085
## Wtr 2 293.147 146.573 18.4591 0.001617 **
## Residuals 7 55.583 7.940
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response tg :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 827.2 206.8 0.6785 0.628246
## Wtr 2 12532.6 6266.3 20.5603 0.001174 **
## Residuals 7 2133.4 304.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response w100 :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 6.974 1.7435 0.3727 0.8214
## Wtr 2 6.200 3.1001 0.6628 0.5450
## Residuals 7 32.741 4.6773
##
## Response Kconc :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 15.3294 3.8323 1.1674 0.4013
## Wtr 2 6.3782 3.1891 0.9715 0.4243
## Residuals 7 22.9792 3.2827
##
## Response pvp :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 0.0019899 0.00049748 1.1446 0.40963
## Wtr 2 0.0030369 0.00151845 3.4936 0.08867 .
## Residuals 7 0.0030425 0.00043464
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Multivariate anova.
Anova(m0)##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 1.6011 0.80092 20 24 0.69052
## Wtr 2 1.5372 2.65685 10 8 0.08971 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(m0, test = "Roy")##
## Type II MANOVA Tests: Roy test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 1.8348 2.2017 5 6 0.182270
## Wtr 2 22.6155 18.0924 5 4 0.007515 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#-----------------------------------------------------------------------
# Residuals.
# Raw residuals.
r <- residuals(m0)
cor(r)## yield tg w100 Kconc pvp
## yield 1.0000000 -0.2721037 0.8476657 -0.16893363 0.24906977
## tg -0.2721037 1.0000000 -0.7317551 0.27715516 0.15051136
## w100 0.8476657 -0.7317551 1.0000000 -0.25482286 0.10770055
## Kconc -0.1689336 0.2771552 -0.2548229 1.00000000 0.08324254
## pvp 0.2490698 0.1505114 0.1077005 0.08324254 1.00000000# Change the 3rd color of the palette used for the ellipses.
oldpal <- palette()
palette(c("black", "blue", "red"))
scatterplotMatrix(r,
gap = 0,
smooth = FALSE,
reg.line = FALSE,
ellipse = TRUE,
diagonal = "qqplot")
palette(oldpal)
corrplot(cor(r),
type = "upper",
tl.pos = "d",
outline = TRUE,
method = "ellipse")
#-----------------------------------------------------------------------
# `wtr` is a 3 level numerical factor.
m1 <- lm(cbind(yield, tg, w100, Kconc, pvp) ~ blk + wtr + I(wtr^2),
data = da)
# Multivariate anova.
anova(m1) # Sequential (more appropriate for the current case).## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99989 5704.7 5 3 3.661e-06 ***
## blk 4 1.49967 0.7 20 24 0.77062
## wtr 1 0.95761 13.6 5 3 0.02851 *
## I(wtr^2) 1 0.58407 0.8 5 3 0.59673
## Residuals 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(m1) # Marginal.##
## Type II MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 1.60110 0.80092 20 24 0.6905
## wtr 1 0.62133 0.98451 5 3 0.5408
## I(wtr^2) 1 0.58407 0.84256 5 3 0.5967anova(m1, test = "Roy") # Sequential (more appropriate ...).## Analysis of Variance Table
##
## Df Roy approx F num Df den Df Pr(>F)
## (Intercept) 1 9507.9 5704.7 5 3 3.661e-06 ***
## blk 4 2.6 3.1 5 6 0.10191
## wtr 1 22.6 13.6 5 3 0.02851 *
## I(wtr^2) 1 1.4 0.8 5 3 0.59673
## Residuals 7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Anova(m1, test = "Roy") # Marginal.##
## Type II MANOVA Tests: Roy test statistic
## Df test stat approx F num Df den Df Pr(>F)
## blk 4 1.8348 2.20172 5 6 0.1823
## wtr 1 1.6409 0.98451 5 3 0.5408
## I(wtr^2) 1 1.4043 0.84256 5 3 0.5967# NOTE: Univariate and Multivariate in disagreement! The way of
# representing the factors (numerical or categorical) is important. BUT
# in this case, the difference is due the type of hypothesis display
# (type I, II or III).
#-----------------------------------------------------------------------
# Performing univariate analysis.
tb <- lapply(c("yield", "tg", "w100", "Kconc", "pvp"),
FUN = function(y) {
m0 <- lm(as.formula(sprintf("%s ~ blk + Wtr", y)),
data = da)
m <- LSmeans(m0,
effect = "Wtr")
g <- glht(model = m0,
linfct = mcp(Wtr = "Tukey"))
g <- cld(g, decreasing = TRUE)$mcletters$Letters
r <- data.frame(
Wtr = m$grid$Wtr,
means = sprintf("%0.2f %s",
m$coef[, "estimate"],
g[as.character(m$grid$Wtr)]))
names(r)[2] <- y
return(r)
})
# Table of means (recursive merging to put the results together).
tbm <- Reduce(merge, x = tb)
tbm## Wtr yield tg w100 Kconc pvp
## 1 37.5 25.31 b 166.80 b 15.24 a 19.61 a 0.99 a
## 2 50 29.79 b 214.40 a 13.99 a 19.01 a 0.96 a
## 3 62.5 37.09 a 241.35 a 15.55 a 17.88 a 1.00 a#-----------------------------------------------------------------------
# Canonical discriminant analysis.
cd_W <- candisc(m0, term = "Wtr")
cd_W##
## Canonical Discriminant Analysis for Wtr:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.95765 22.6155 21.237 94.2564 94.256
## 2 0.57950 1.3781 21.237 5.7436 100.000
##
## Test of H0: The canonical correlations in the
## current row and all that follow are zero
##
## LR test stat approx F numDF denDF Pr(> F)
## 1 0.01781 9.0916 10 14 0.0001516 ***
## 2 0.42050 2.7562 4 8 0.1037425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(cd_W)##
## Canonical Discriminant Analysis for Wtr:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.95765 22.6155 21.237 94.2564 94.256
## 2 0.57950 1.3781 21.237 5.7436 100.000
##
## Class means:
##
## Can1 Can2
## 37.5 4.04538 0.52335
## 50 -0.50224 -1.08863
## 62.5 -4.42892 0.70660
##
## std coefficients:
## Can1 Can2
## yield -4.13965 2.405074
## tg 1.88860 -2.261965
## w100 4.89529 -3.432387
## Kconc 0.11646 -0.020264
## pvp 0.19559 0.961866# The weights to compute the scores.
cd_W$coeffs.raw## Can1 Can2
## yield -1.46906820 0.85350564
## tg 0.10818053 -0.12956727
## w100 2.26350008 -1.58707917
## Kconc 0.06427546 -0.01118409
## pvp 9.38174341 46.13703065# The scores.
head(cd_W$scores)## blk Wtr Can1 Can2
## 4 I 37.5 3.750538 0.1245453
## 9 I 50 -1.988107 -2.2199431
## 14 I 62.5 -6.303187 0.4007003
## 19 II 37.5 5.794425 0.6828719
## 24 II 50 1.003528 0.6133912
## 29 II 62.5 -4.684761 0.6964649# Calculating the scores (with the centered y-variables).
Y <- m0$model[[1]]
z <- scale(Y, scale = FALSE) %*% cd_W$coeffs.raw
head(z[, 1:2])## Can1 Can2
## [1,] 3.750538 0.1245453
## [2,] -1.988107 -2.2199431
## [3,] -6.303187 0.4007003
## [4,] 5.794425 0.6828719
## [5,] 1.003528 0.6133912
## [6,] -4.684761 0.6964649head(cd_W$scores[, c("Can1", "Can2")])## Can1 Can2
## 4 3.750538 0.1245453
## 9 -1.988107 -2.2199431
## 14 -6.303187 0.4007003
## 19 5.794425 0.6828719
## 24 1.003528 0.6133912
## 29 -4.684761 0.6964649# The loadings are the correlations between scores and responses.
cd_W$structure## Can1 Can2
## yield -0.89843778 0.17624288
## tg -0.93341861 -0.20325992
## w100 0.04883054 0.50377294
## Kconc 0.38183868 0.09871581
## pvp 0.06725178 0.91739356cor(cbind(Y, z))[colnames(Y), colnames(z)]## Can1 Can2
## yield -0.89843778 0.17624288
## tg -0.93341861 -0.20325992
## w100 0.04883054 0.50377294
## Kconc 0.38183868 0.09871581
## pvp 0.06725178 0.91739356# The biplot.
plot(cd_W)## Vector scale factor set to 4.946
2.3 Analysing the entire soybeanwp dataset
#----------------------------------------------------------------------
# Fitting the multivariate Gaussian linear model.
xyplot(yield + tg + w100 + Kconc + pvp ~ pts,
groups = Wtr,
auto.key = TRUE,
type = c("p", "a"),
outer = TRUE,
as.table = TRUE,
data = swp[-74, ],
xlab = "Potassium",
ylab = "Values",
scales = list(y = "free"))
m0 <- lm(cbind(yield, tg, w100, Kconc, pvp) ~ blk + Wtr * Pts,
data = swp[-74, ])
# MANOVA.
anova(m0)## Analysis of Variance Table
##
## Df Pillai approx F num Df den Df Pr(>F)
## (Intercept) 1 0.99979 47622 5 51 < 2.2e-16 ***
## blk 4 0.55522 2 20 216 0.0289970 *
## Wtr 2 1.04187 11 10 104 6.971e-13 ***
## Pts 4 1.54050 7 20 216 2.732e-14 ***
## Wtr:Pts 8 1.17190 2 40 275 0.0002756 ***
## Residuals 55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# ANOVAs.
summary.aov(m0)## Response yield :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 100.39 25.10 4.8948 0.001896 **
## Wtr 2 446.53 223.26 43.5456 4.617e-12 ***
## Pts 4 2592.13 648.03 126.3925 < 2.2e-16 ***
## Wtr:Pts 8 286.00 35.75 6.9726 2.629e-06 ***
## Residuals 55 281.99 5.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response tg :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 5108 1277.1 2.9635 0.027419 *
## Wtr 2 32505 16252.6 37.7139 4.867e-11 ***
## Pts 4 65223 16305.8 37.8374 3.421e-15 ***
## Wtr:Pts 8 11397 1424.6 3.3059 0.003748 **
## Residuals 55 23702 430.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response w100 :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 8.019 2.0046 1.0335 0.398303
## Wtr 2 26.220 13.1099 6.7588 0.002374 **
## Pts 4 122.229 30.5573 15.7537 1.194e-08 ***
## Wtr:Pts 8 33.609 4.2011 2.1659 0.044523 *
## Residuals 55 106.683 1.9397
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response Kconc :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 13.59 3.399 0.9782 0.4270
## Wtr 2 7.03 3.517 1.0122 0.3701
## Pts 4 448.68 112.171 32.2865 7.503e-14 ***
## Wtr:Pts 8 17.13 2.141 0.6162 0.7605
## Residuals 55 191.08 3.474
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Response pvp :
## Df Sum Sq Mean Sq F value Pr(>F)
## blk 4 0.00391 0.000977 0.6282 0.644399
## Wtr 2 0.01937 0.009687 6.2286 0.003645 **
## Pts 4 0.48955 0.122388 78.6930 < 2.2e-16 ***
## Wtr:Pts 8 0.02342 0.002928 1.8823 0.081480 .
## Residuals 55 0.08554 0.001555
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Univariate model summaries.
# summary(m0)
# Estimated parameter.
coef(m0)## yield tg w100 Kconc pvp
## (Intercept) 14.0727143 125.088095 11.2442738 14.3074048 0.796763719
## blkII 1.0660000 -3.666667 0.9466667 0.8626667 0.022717159
## blkIII -0.8606667 -6.866667 0.1793333 0.6960000 0.011758773
## blkIV -1.4433333 -19.133333 0.5713333 0.6640000 0.009871475
## blkV -1.5255714 -16.773810 0.5412976 1.3903095 0.015282874
## Wtr50 2.1920000 16.400000 0.1940000 -1.7960000 -0.082674174
## Wtr62.5 2.5700000 23.600000 -0.1480000 -2.2900000 -0.045067147
## Pts30 6.8140000 40.200000 1.3340000 1.9900000 0.142767113
## Pts60 10.4060000 47.600000 3.1680000 2.6880000 0.175568777
## Pts120 11.7880000 51.000000 3.5480000 4.5820000 0.180467734
## Pts180 11.8700000 39.800000 4.7500000 5.8740000 0.180873826
## Wtr50:Pts30 0.0440000 13.600000 -1.0860000 1.7960000 0.070467482
## Wtr62.5:Pts30 -1.9160000 -15.800000 -0.0220000 1.7920000 -0.027963182
## Wtr50:Pts60 2.5740000 29.600000 -1.2580000 3.0900000 0.061348687
## Wtr62.5:Pts60 3.3320000 28.600000 -0.8800000 2.8880000 0.027629784
## Wtr50:Pts120 2.2860000 31.200000 -1.4420000 1.1940000 0.055494778
## Wtr62.5:Pts120 8.4612857 46.478571 0.3338929 1.0849286 0.048622758
## Wtr50:Pts180 1.1780000 65.000000 -4.1820000 2.0920000 0.040734774
## Wtr62.5:Pts180 9.1860000 53.200000 -0.2280000 2.2900000 0.042618230# Raw residuals.
r <- residuals(m0)
var(r)## yield tg w100 Kconc pvp
## yield 3.86291731 17.0772655 8.817030e-01 -1.297742387 -0.0076469699
## tg 17.07726549 324.6840835 -1.335361e+01 -7.802964775 -0.0639091046
## w100 0.88170299 -13.3536070 1.461406e+00 -0.206366513 -0.0007597463
## Kconc -1.29774239 -7.8029648 -2.063665e-01 2.617582172 0.0043406883
## pvp -0.00764697 -0.0639091 -7.597463e-04 0.004340688 0.0011717687cor(r)## yield tg w100 Kconc pvp
## yield 1.0000000 0.4822032 0.37108989 -0.40811318 -0.11366080
## tg 0.4822032 1.0000000 -0.61303124 -0.26765723 -0.10361229
## w100 0.3710899 -0.6130312 1.00000000 -0.10551242 -0.01835956
## Kconc -0.4081132 -0.2676572 -0.10551242 1.00000000 0.07837679
## pvp -0.1136608 -0.1036123 -0.01835956 0.07837679 1.00000000# Change the 3rd color of the palette used for the ellipses.
oldpal <- palette()
palette(c("black", "blue", "red"))
scatterplotMatrix(r,
gap = 0,
smooth = FALSE,
reg.line = FALSE,
ellipse = TRUE,
diagonal = "qqplot")
palette(oldpal)
corrplot(cor(r),
type = "upper",
tl.pos = "d",
outline = TRUE,
method = "ellipse")
#-----------------------------------------------------------------------
# Canonical analysis.
cd_WP <- candisc(m0, term = "Wtr:Pts")
cd_WP##
## Canonical Discriminant Analysis for Wtr:Pts:
##
## CanRsq Eigenvalue Difference Percent Cumulative
## 1 0.690657 2.232655 1.8329 78.42318 78.423
## 2 0.285574 0.399725 1.8329 14.04053 92.464
## 3 0.112497 0.126757 1.8329 4.45241 96.916
## 4 0.062688 0.066880 1.8329 2.34921 99.265
## 5 0.020487 0.020916 1.8329 0.73467 100.000
##
## Test of H0: The canonical correlations in the
## current row and all that follow are zero
##
## LR test stat approx F numDF denDF Pr(> F)
## 1 0.18008 3.2368 40 268.69 7.348e-09 ***
## 2 0.58213 1.3008 28 224.97 0.1513
## 3 0.81482 0.7453 18 178.68 0.7606
## 4 0.91811 0.5587 10 128.00 0.8448
## 5 0.97951 0.3399 4 65.00 0.8501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Correlations between responses and canonical scores.
cd_WP$structure## Can1 Can2 Can3 Can4 Can5
## yield -0.9675764 -0.6369567 0.6255565 0.02493008 -0.8279715
## tg -0.8259545 -0.8312523 0.4255786 -0.15046779 -0.6774879
## w100 -0.5550029 -0.0481521 0.6673181 0.34247061 -0.6586667
## Kconc -0.6274747 -0.6606093 0.5946264 0.75127527 -0.8309136
## pvp -0.7246166 -0.7093252 0.9450160 0.59791723 -0.7352603# Canonical scores.
head(cd_WP$scores)## blk Wtr Pts Can1 Can2 Can3 Can4 Can5
## 1 I 37.5 0 4.5062699 2.64196284 -5.4440764 -0.9498030 4.220891
## 2 I 37.5 30 1.2441462 0.08125175 1.1742776 0.6982192 1.105535
## 3 I 37.5 60 0.3028206 0.06520864 2.5073658 1.2613717 -1.376520
## 4 I 37.5 120 0.7034910 -0.55018488 0.8036896 0.9272220 0.910372
## 5 I 37.5 180 -1.9652384 -0.58786511 1.6976880 0.7285266 -1.090994
## 6 I 50 0 4.6560107 2.72948577 -3.1448520 -2.1000088 3.821640# Biplot.
plot(cd_WP)## Vector scale factor set to 9.221
# Merge the responses with the scores.
swps <- merge(swp, cd_WP$scores)
xyplot(-Can2 ~ -Can1 | Wtr,
groups = Pts,
auto.key = TRUE,
data = swps) +
glayer(panel.ellipse(...))
xyplot(-Can1 + -Can2 + yield + tg + w100 + Kconc ~ Pts,
groups = Wtr,
outer = TRUE,
scales = list(y = "free"),
as.table = TRUE,
data = swps,
type = c("p", "a"))
3 Multivariate analysis using the mcglm package
- Five outcomes yield, w100, Kconc (continuous), tg (count) and nvp (binomial).
- Modelling counting outcome (tg).
# Linear predictor
f.tg <- tg ~ block + water * potassium
# Matrix linear predictor
Z0 <- mc_id(soybeanwp)
# Fitting extended Poisson-Tweedie models
fit.tg <- mcglm(linear_pred = c(f.tg), matrix_pred = list(c(Z0)),
link = "log", variance = "poisson_tweedie",
power_fixed = FALSE,
data = soybeanwp)## Automatic initial values selected.summary(fit.tg, print = c("power","Dispersion"))## Power:
## Estimates Std.error Z value
## 1 0.5382939 0.8852954 0.6080388
##
## Dispersion:
## Estimates Std.error Z value
## 1 36.90082 168.4576 0.219051
##
## Algorithm: chaser
## Correction: TRUE
## Number iterations: 19anova(fit.tg)## Wald test for fixed effects
## Call: tg ~ block + water * potassium
##
## Covariate Chi.Square Df p.value
## 1 blockII 8.2322 4 0.0834
## 2 water 2.8491 1 0.0914
## 3 potassium 0.1112 1 0.7387
## 4 water:potassium 3.3405 1 0.0676- Modelling the binomial outcome.
# Computing the proportions
Ntrial <- soybeanwp$nip + soybeanwp$nvp
soybeanwp$prop <- soybeanwp$nvp/(Ntrial)
# Linear predictor
f.nvp <- prop ~ block + water * potassium
# Matrix linear predictor
Z0 <- mc_id(soybeanwp)
# Fitting extended Poisson-Tweedie models
fit.nvp <- mcglm(linear_pred = c(f.nvp), matrix_pred = list(c(Z0)),
link = "logit", variance = "binomialP",
power_fixed = TRUE, Ntrial = list(Ntrial),
data = soybeanwp)## Automatic initial values selected.summary(fit.nvp, print = c("Dispersion"))## Dispersion:
## Estimates Std.error Z value
## 1 3.743841 1.189834 3.146522
##
## Algorithm: chaser
## Correction: TRUE
## Number iterations: 7anova(fit.nvp)## Wald test for fixed effects
## Call: prop ~ block + water * potassium
##
## Covariate Chi.Square Df p.value
## 1 blockII 1.5759 4 0.8131
## 2 water 1.4559 1 0.2276
## 3 potassium 1.0955 1 0.2953
## 4 water:potassium 0.1845 1 0.6675- Modelling the continuous outcomes.
# Linear predictors
f.yield <- yield ~ block + water * potassium
f.w100 <- w100 ~ block + water * potassium
f.Kconc <- Kconc ~ block + water * potassium
# Fitting univariate Gaussian models
# Outcome: yield
fit1.yield <- mcglm(linear_pred = c(f.yield), matrix_pred = list(c(Z0)),
data = soybeanwp)## Automatic initial values selected.# Outcome: w100
fit1.w100 <- mcglm(linear_pred = c(f.w100), matrix_pred = list(c(Z0)),
data = soybeanwp)## Automatic initial values selected.# Outcome: Kconc
fit1.Kconc <- mcglm(linear_pred = c(f.Kconc), matrix_pred = list(c(Z0)),
data = soybeanwp)## Automatic initial values selected.# MANOVA for the three continuous outcomes
manova_fit <- manova(cbind(yield, w100, Kconc) ~ block + water * potassium,
data = soybeanwp)
# MANOVA fitted by the mcglm package
manova_mcglm <- mcglm(linear_pred = c(f.yield, f.w100, f.Kconc),
matrix_pred = list(c(Z0),c(Z0),c(Z0)),
control_algorithm = list(correct = FALSE),
data = soybeanwp)## Automatic initial values selected.# Comparing the fits (they provide exactly the same fit)
source("../review/functions.R")
logLik.mlm(manova_fit) # MANOVA function## 'log Lik.' -494.1584 (df=30)plogLik(manova_mcglm) # McGLM function## Pseudo log Lik. -494.16 (df=30)# Multivariate hypotheses tests
Anova(manova_fit, test = "Hotelling-Lawley", type = "III")##
## Type III MANOVA Tests: Hotelling-Lawley test statistic
## Df test stat approx F num Df den Df Pr(>F)
## (Intercept) 1 2.11836 45.898 3 65 4.785e-16 ***
## block 4 0.24531 1.302 12 191 0.220269
## water 1 0.05020 1.088 3 65 0.360697
## potassium 1 0.08479 1.837 3 65 0.149166
## water:potassium 1 0.19356 4.194 3 65 0.008937 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1manova.mcglm(manova_mcglm)## Effects Df Hotelling.Lawley Qui.square p_value
## 1 Intercept 3 2.118 158.877 3.202926e-34
## 2 block 12 0.245 18.398 1.041243e-01
## 3 water 3 0.050 3.765 2.880228e-01
## 4 potassium 3 0.085 6.359 9.539460e-02
## 5 water:potassium 3 0.194 14.517 2.279408e-03- What about the counting and binomial outcomes?
mcglm_4all <- mcglm(linear_pred = c(f.yield, f.w100, f.Kconc, f.tg, f.nvp),
matrix_pred = list(c(Z0),c(Z0),c(Z0),c(Z0), c(Z0)),
link = c("identity","identity","identity","log","logit"),
variance = c("constant","constant","constant",
"poisson_tweedie","binomialP"),
Ntrial = list(NULL, NULL, NULL, NULL, Ntrial),
control_algorithm = list(correct = FALSE),
data = soybeanwp)## Automatic initial values selected.# Multivariate hypotheses tests for non-Gaussian data
manova.mcglm(mcglm_4all)## Effects Df Hotelling.Lawley Qui.square p_value
## 1 Intercept 5 280.387 21029.043 0.000000e+00
## 2 block 20 0.248 18.610 5.472706e-01
## 3 water 5 0.152 11.406 4.390387e-02
## 4 potassium 5 0.674 50.551 1.068983e-09
## 5 water:potassium 5 0.803 60.200 1.104905e-11# Correlation matrix between responses
COR <- Matrix(0, 5, 5)
COR[lower.tri(COR)] <- mcglm_4all$Covariance[1:10]
COR <- matrix(forceSymmetric(COR, uplo = TRUE), 5, 5)
colnames(COR) <- rownames(COR) <- c("yield", "w100", "K", "tg", "nvp")
corrplot(COR, type = "upper",
tl.pos = "d",
outline = TRUE,
method = "ellipse")
4 References
Serafim, M.E., Ono, F.B., Zeviani, W.M., Novelino, J.O., Silva, J.V., 2012. Umidade do solo e doses de potássio na cultura da soja. Revista Ciência Agronômica 43, 222–227.
5 Session information
# devtools::session_info()
Sys.time()## [1] "2017-07-26 19:36:08 BRT"cbind(Sys.info())## [,1]
## sysname "Linux"
## release "4.4.0-87-generic"
## version "#110-Ubuntu SMP Tue Jul 18 12:55:35 UTC 2017"
## nodename "class"
## machine "x86_64"
## login "unknown"
## user "walmes"
## effective_user "walmes"sessionInfo()## R version 3.4.1 (2017-06-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
##
## 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] methods stats graphics grDevices utils datasets base
##
## other attached packages:
## [1] Matrix_1.2-10 mcglm_0.4.0 multcomp_1.4-5
## [4] TH.data_1.0-8 MASS_7.3-47 survival_2.40-1
## [7] mvtnorm_1.0-5 doBy_4.5-15 corrplot_0.77
## [10] candisc_0.7-2 heplots_1.3-3 car_2.1-4
## [13] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-35
## [16] knitr_1.15.1 rmarkdown_1.3
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.9 compiler_3.4.1 nloptr_1.0.4
## [4] tools_3.4.1 digest_0.6.12 lme4_1.1-12
## [7] evaluate_0.10 nlme_3.1-131 mgcv_1.8-17
## [10] yaml_2.1.14 parallel_3.4.1 SparseM_1.74
## [13] stringr_1.2.0 MatrixModels_0.4-1 rprojroot_1.2
## [16] grid_3.4.1 nnet_7.3-12 minqa_1.2.4
## [19] magrittr_1.5 codetools_0.2-15 backports_1.0.5
## [22] htmltools_0.3.5 splines_3.4.1 assertthat_0.2.0
## [25] pbkrtest_0.4-6 quantreg_5.29 sandwich_2.3-4
## [28] stringi_1.1.2 zoo_1.7-14
mcglm4aed - Multivariate Covariance Generalized
Linear Models for the Analysis of Experimental Data
|
Bonat & Zeviani (2017) |