|
Extensões de Modelos de Regressão
|
Tópicos em Modelos Mistos
Prof. Walmes M. Zeviani & Prof. Paulo Justiniano R. Jr
2017-09-04
1 Delineamento de blocos completos casualizados
# Lista de pacotes.
pkg <- c("lattice",
"latticeExtra",
"gridExtra",
"reshape",
"car",
"lme4",
"doBy",
"multcomp",
"contrast",
"wzRfun",
"labestData")
sapply(pkg,
FUN = library,
character.only = TRUE,
logical.return = TRUE)1.1 Modelo de efeitos fixos
# da <- read.table("http://www.leg.ufpr.br/~walmes/data/pimentel_batatinha.txt",
# header = TRUE, sep = "\t")
# labestDataView()
# help(PimentelEg5.2, help_type = "html")
xtabs(~variedade + bloco, data = PimentelEg5.2)## bloco
## variedade I II III IV
## B 116-51 1 1 1 1
## B 1-52 1 1 1 1
## B 25-50 E 1 1 1 1
## B 72-53 A 1 1 1 1
## Buena Vista 1 1 1 1
## Huinkul 1 1 1 1
## Kennebec 1 1 1 1
## S. Rafalela 1 1 1 1xyplot(producao ~ variedade,
groups = bloco,
data = PimentelEg5.2,
type = "o",
ylab = expression(Produção ~ (t ~ ha^{-1})),
xlab = "Variedades de batatinha")
da <- PimentelEg5.2
str(da)## 'data.frame': 32 obs. of 3 variables:
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 1 1 1 1 1 1 1 2 2 ...
## $ variedade: Factor w/ 8 levels "B 116-51","B 1-52",..: 1 2 3 4 5 6 7 8 1 2 ...
## $ producao : num 23.1 20 12.7 18 15.4 21.1 9.2 22.6 24.2 21.1 ...#-----------------------------------------------------------------------
# Ajuste o modelo linear de efeitos fixos.
m0 <- lm(producao ~ bloco + variedade, data = da)
# Diagnóstico.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: producao
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 50.53 16.843 1.9709 0.1493
## variedade 7 919.72 131.389 15.3744 5.723e-07 ***
## Residuals 21 179.47 8.546
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Medidas de ajuste e estimativas dos parâmetros (efeitos).
summary(m0)##
## Call:
## lm(formula = producao ~ bloco + variedade, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2750 -1.6437 0.0625 1.0563 5.6500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.550 1.714 11.990 7.40e-11 ***
## blocoII 3.500 1.462 2.395 0.02605 *
## blocoIII 2.275 1.462 1.556 0.13455
## blocoIV 2.025 1.462 1.385 0.18047
## variedadeB 1-52 -0.225 2.067 -0.109 0.91436
## variedadeB 25-50 E -6.000 2.067 -2.903 0.00851 **
## variedadeB 72-53 A 0.300 2.067 0.145 0.88599
## variedadeBuena Vista -10.075 2.067 -4.874 8.08e-05 ***
## variedadeHuinkul 2.550 2.067 1.234 0.23098
## variedadeKennebec -11.800 2.067 -5.708 1.15e-05 ***
## variedadeS. Rafalela 2.950 2.067 1.427 0.16825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.923 on 21 degrees of freedom
## Multiple R-squared: 0.8439, Adjusted R-squared: 0.7696
## F-statistic: 11.35 on 10 and 21 DF, p-value: 2.153e-06#-----------------------------------------------------------------------
# Médias ajustadas, usando contrast::contrast().
# help(contrast, help_type = "html")
# Primeira variedade no primeiro bloco.
contrast(fit = m0,
a = list(variedade = levels(m0$model$variedade)[1],
bloco = levels(m0$model$bloco)[1]))## lm model parameter contrast
##
## Contrast S.E. Lower Upper t df Pr(>|t|)
## 20.55 1.713964 16.98562 24.11438 11.99 21 0# Primeira variedade na média do efeito de bloco.
contrast(fit = m0,
a = list(variedade = levels(m0$model$variedade)[1],
bloco = levels(m0$model$bloco)),
type = "average")## lm model parameter contrast
##
## Contrast S.E. Lower Upper t df Pr(>|t|)
## 1 22.5 1.461673 19.46028 25.53972 15.39 21 0# Todas as variedades no primeiro bloco.
contrast(fit = m0,
a = list(variedade = levels(m0$model$variedade),
bloco = levels(m0$model$bloco)[1]))## lm model parameter contrast
##
## Contrast S.E. Lower Upper t df Pr(>|t|)
## 20.550 1.713964 16.985618 24.11438 11.99 21 0
## 20.325 1.713964 16.760618 23.88938 11.86 21 0
## 14.550 1.713964 10.985618 18.11438 8.49 21 0
## 20.850 1.713964 17.285618 24.41438 12.16 21 0
## 10.475 1.713964 6.910618 14.03938 6.11 21 0
## 23.100 1.713964 19.535618 26.66438 13.48 21 0
## 8.750 1.713964 5.185618 12.31438 5.11 21 0
## 23.500 1.713964 19.935618 27.06438 13.71 21 0#-----------------------------------------------------------------------
# Estudando a estrutura do objeto.
ctr <- contrast(fit = m0,
a = list(variedade = levels(m0$model$variedade),
bloco = levels(m0$model$bloco)[1]))
names(ctr)## [1] "variedade" "bloco" "Contrast" "SE" "Lower"
## [6] "Upper" "testStat" "df" "Pvalue" "var"
## [11] "X" "df.residual" "nvary" "foldChange" "aCoef"
## [16] "bCoef" "model" "covType"t(ctr$X)## 1 2 3 4 5 6 7 8
## (Intercept) 1 1 1 1 1 1 1 1
## blocoII 0 0 0 0 0 0 0 0
## blocoIII 0 0 0 0 0 0 0 0
## blocoIV 0 0 0 0 0 0 0 0
## variedadeB 1-52 0 1 0 0 0 0 0 0
## variedadeB 25-50 E 0 0 1 0 0 0 0 0
## variedadeB 72-53 A 0 0 0 1 0 0 0 0
## variedadeBuena Vista 0 0 0 0 1 0 0 0
## variedadeHuinkul 0 0 0 0 0 1 0 0
## variedadeKennebec 0 0 0 0 0 0 1 0
## variedadeS. Rafalela 0 0 0 0 0 0 0 1
## attr(,"assign")
## [1] 0 1 1 1 2 2 2 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$bloco
## [1] "contr.treatment"
##
## attr(,"contrasts")$variedade
## [1] "contr.treatment"# Ponderar nas entradas dos pesos de blocos.
ctr$X[, m0$assign == 1] <- 1/nlevels(da$bloco)
t(ctr$X)## 1 2 3 4 5 6 7 8
## (Intercept) 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
## blocoII 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
## blocoIII 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
## blocoIV 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
## variedadeB 1-52 0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00
## variedadeB 25-50 E 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00
## variedadeB 72-53 A 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00
## variedadeBuena Vista 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00
## variedadeHuinkul 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00
## variedadeKennebec 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00
## variedadeS. Rafalela 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
## attr(,"assign")
## [1] 0 1 1 1 2 2 2 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$bloco
## [1] "contr.treatment"
##
## attr(,"contrasts")$variedade
## [1] "contr.treatment"# Médias ajustadas.
ctr$X %*% coef(m0)## [,1]
## 1 22.500
## 2 22.275
## 3 16.500
## 4 22.800
## 5 12.425
## 6 25.050
## 7 10.700
## 8 25.450# Matriz de covariância das médias.
round(ctr$X %*% vcov(m0) %*% t(ctr$X), digits = 4)## 1 2 3 4 5 6 7 8
## 1 2.1365 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## 2 0.0000 2.1365 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
## 3 0.0000 0.0000 2.1365 0.0000 0.0000 0.0000 0.0000 0.0000
## 4 0.0000 0.0000 0.0000 2.1365 0.0000 0.0000 0.0000 0.0000
## 5 0.0000 0.0000 0.0000 0.0000 2.1365 0.0000 0.0000 0.0000
## 6 0.0000 0.0000 0.0000 0.0000 0.0000 2.1365 0.0000 0.0000
## 7 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2.1365 0.0000
## 8 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2.1365#-----------------------------------------------------------------------
# Usando o doBy::LSmeans().
# help(LSmeans, help_type = "html")
# Médias das variedades considerando o efeito médio dos blocos.
LSmeans(object = m0, effect = "variedade")## estimate se df t.stat p.value variedade
## 1 22.500 1.461673 21 15.393319 6.524482e-13 B 116-51
## 2 22.275 1.461673 21 15.239386 7.925028e-13 B 1-52
## 3 16.500 1.461673 21 11.288434 2.228887e-10 B 25-50 E
## 4 22.800 1.461673 21 15.598564 5.047112e-13 B 72-53 A
## 5 12.425 1.461673 21 8.500533 3.070928e-08 Buena Vista
## 6 25.050 1.461673 21 17.137896 8.027793e-14 Huinkul
## 7 10.700 1.461673 21 7.320379 3.315981e-07 Kennebec
## 8 25.450 1.461673 21 17.411555 5.877653e-14 S. Rafalela# Fixando no efeito do primeiro bloco.
LSmeans(object = m0, effect = "variedade",
at = list(bloco = levels(da$bloco)[1]))## estimate se df t.stat p.value variedade bloco
## 1 20.550 1.713964 21 11.989753 7.397251e-11 B 116-51 I
## 2 20.325 1.713964 21 11.858478 9.059641e-11 B 1-52 I
## 3 14.550 1.713964 21 8.489095 3.139780e-08 B 25-50 E I
## 4 20.850 1.713964 21 12.164786 5.659946e-11 B 72-53 A I
## 5 10.475 1.713964 21 6.111565 4.594416e-06 Buena Vista I
## 6 23.100 1.713964 21 13.477533 8.314969e-12 Huinkul I
## 7 8.750 1.713964 21 5.105126 4.678097e-05 Kennebec I
## 8 23.500 1.713964 21 13.710910 6.005879e-12 S. Rafalela Ilsm <- LSmeans(object = m0, effect = "variedade")
t(lsm$K)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## (Intercept) 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
## blocoII 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
## blocoIII 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
## blocoIV 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
## variedadeB 1-52 0.00 1.00 0.00 0.00 0.00 0.00 0.00 0.00
## variedadeB 25-50 E 0.00 0.00 1.00 0.00 0.00 0.00 0.00 0.00
## variedadeB 72-53 A 0.00 0.00 0.00 1.00 0.00 0.00 0.00 0.00
## variedadeBuena Vista 0.00 0.00 0.00 0.00 1.00 0.00 0.00 0.00
## variedadeHuinkul 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00
## variedadeKennebec 0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00
## variedadeS. Rafalela 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00#-----------------------------------------------------------------------
# Comparações múltiplas de médias.
# A função cria a matriz de contrastes fixando os fatores categóricos
# não mencionados no nível de referência, no caso, bloco I.
g0 <- summary(glht(m0, linfct = mcp(variedade = "Tukey")),
test = adjusted(type = "fdr"))
g0##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = producao ~ bloco + variedade, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B 1-52 - B 116-51 == 0 -0.225 2.067 -0.109 0.914357
## B 25-50 E - B 116-51 == 0 -6.000 2.067 -2.903 0.017029 *
## B 72-53 A - B 116-51 == 0 0.300 2.067 0.145 0.914357
## Buena Vista - B 116-51 == 0 -10.075 2.067 -4.874 0.000251 ***
## Huinkul - B 116-51 == 0 2.550 2.067 1.234 0.293976
## Kennebec - B 116-51 == 0 -11.800 2.067 -5.708 5.36e-05 ***
## S. Rafalela - B 116-51 == 0 2.950 2.067 1.427 0.247947
## B 25-50 E - B 1-52 == 0 -5.775 2.067 -2.794 0.019042 *
## B 72-53 A - B 1-52 == 0 0.525 2.067 0.254 0.898222
## Buena Vista - B 1-52 == 0 -9.850 2.067 -4.765 0.000293 ***
## Huinkul - B 1-52 == 0 2.775 2.067 1.342 0.271297
## Kennebec - B 1-52 == 0 -11.575 2.067 -5.600 5.91e-05 ***
## S. Rafalela - B 1-52 == 0 3.175 2.067 1.536 0.216969
## B 72-53 A - B 25-50 E == 0 6.300 2.067 3.048 0.013172 *
## Buena Vista - B 25-50 E == 0 -4.075 2.067 -1.971 0.102125
## Huinkul - B 25-50 E == 0 8.550 2.067 4.136 0.001095 **
## Kennebec - B 25-50 E == 0 -5.800 2.067 -2.806 0.019042 *
## S. Rafalela - B 25-50 E == 0 8.950 2.067 4.330 0.000752 ***
## Buena Vista - B 72-53 A == 0 -10.375 2.067 -5.019 0.000201 ***
## Huinkul - B 72-53 A == 0 2.250 2.067 1.088 0.351487
## Kennebec - B 72-53 A == 0 -12.100 2.067 -5.854 4.62e-05 ***
## S. Rafalela - B 72-53 A == 0 2.650 2.067 1.282 0.285098
## Huinkul - Buena Vista == 0 12.625 2.067 6.108 3.25e-05 ***
## Kennebec - Buena Vista == 0 -1.725 2.067 -0.834 0.482294
## S. Rafalela - Buena Vista == 0 13.025 2.067 6.301 2.81e-05 ***
## Kennebec - Huinkul == 0 -14.350 2.067 -6.942 1.04e-05 ***
## S. Rafalela - Huinkul == 0 0.400 2.067 0.194 0.913685
## S. Rafalela - Kennebec == 0 14.750 2.067 7.136 1.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)# A função que passa os contrastes é essa.
contrMat(n = 1:nlevels(da$variedade), type = "Tukey")##
## Multiple Comparisons of Means: Tukey Contrasts
##
## 1 2 3 4 5 6 7 8
## 2 - 1 -1 1 0 0 0 0 0 0
## 3 - 1 -1 0 1 0 0 0 0 0
## 4 - 1 -1 0 0 1 0 0 0 0
## 5 - 1 -1 0 0 0 1 0 0 0
## 6 - 1 -1 0 0 0 0 1 0 0
## 7 - 1 -1 0 0 0 0 0 1 0
## 8 - 1 -1 0 0 0 0 0 0 1
## 3 - 2 0 -1 1 0 0 0 0 0
## 4 - 2 0 -1 0 1 0 0 0 0
## 5 - 2 0 -1 0 0 1 0 0 0
## 6 - 2 0 -1 0 0 0 1 0 0
## 7 - 2 0 -1 0 0 0 0 1 0
## 8 - 2 0 -1 0 0 0 0 0 1
## 4 - 3 0 0 -1 1 0 0 0 0
## 5 - 3 0 0 -1 0 1 0 0 0
## 6 - 3 0 0 -1 0 0 1 0 0
## 7 - 3 0 0 -1 0 0 0 1 0
## 8 - 3 0 0 -1 0 0 0 0 1
## 5 - 4 0 0 0 -1 1 0 0 0
## 6 - 4 0 0 0 -1 0 1 0 0
## 7 - 4 0 0 0 -1 0 0 1 0
## 8 - 4 0 0 0 -1 0 0 0 1
## 6 - 5 0 0 0 0 -1 1 0 0
## 7 - 5 0 0 0 0 -1 0 1 0
## 8 - 5 0 0 0 0 -1 0 0 1
## 7 - 6 0 0 0 0 0 -1 1 0
## 8 - 6 0 0 0 0 0 -1 0 1
## 8 - 7 0 0 0 0 0 0 -1 1# Tipos de contraste disponíveis.
formals(fun = contrMat)## $n
##
##
## $type
## c("Dunnett", "Tukey", "Sequen", "AVE", "Changepoint", "Williams",
## "Marcus", "McDermott", "UmbrellaWilliams", "GrandMean")
##
## $base
## [1] 1# Resumo compacto com letras.
cld(g0)## B 116-51 B 1-52 B 25-50 E B 72-53 A Buena Vista Huinkul
## "c" "c" "b" "c" "ab" "c"
## Kennebec S. Rafalela
## "a" "c"# Criando a matriz de contrastes de Tukey (contrastes par a par).
allctr <- wzRfun::apc(lfm = lsm$K)
allctr## (Intercept) blocoII blocoIII blocoIV
## B 116-51-B 1-52 0 0 0 0
## B 116-51-B 25-50 E 0 0 0 0
## B 116-51-B 72-53 A 0 0 0 0
## B 116-51-Buena Vista 0 0 0 0
## B 116-51-Huinkul 0 0 0 0
## B 116-51-Kennebec 0 0 0 0
## B 116-51-S. Rafalela 0 0 0 0
## B 1-52-B 25-50 E 0 0 0 0
## B 1-52-B 72-53 A 0 0 0 0
## B 1-52-Buena Vista 0 0 0 0
## B 1-52-Huinkul 0 0 0 0
## B 1-52-Kennebec 0 0 0 0
## B 1-52-S. Rafalela 0 0 0 0
## B 25-50 E-B 72-53 A 0 0 0 0
## B 25-50 E-Buena Vista 0 0 0 0
## B 25-50 E-Huinkul 0 0 0 0
## B 25-50 E-Kennebec 0 0 0 0
## B 25-50 E-S. Rafalela 0 0 0 0
## B 72-53 A-Buena Vista 0 0 0 0
## B 72-53 A-Huinkul 0 0 0 0
## B 72-53 A-Kennebec 0 0 0 0
## B 72-53 A-S. Rafalela 0 0 0 0
## Buena Vista-Huinkul 0 0 0 0
## Buena Vista-Kennebec 0 0 0 0
## Buena Vista-S. Rafalela 0 0 0 0
## Huinkul-Kennebec 0 0 0 0
## Huinkul-S. Rafalela 0 0 0 0
## Kennebec-S. Rafalela 0 0 0 0
## variedadeB 1-52 variedadeB 25-50 E
## B 116-51-B 1-52 -1 0
## B 116-51-B 25-50 E 0 -1
## B 116-51-B 72-53 A 0 0
## B 116-51-Buena Vista 0 0
## B 116-51-Huinkul 0 0
## B 116-51-Kennebec 0 0
## B 116-51-S. Rafalela 0 0
## B 1-52-B 25-50 E 1 -1
## B 1-52-B 72-53 A 1 0
## B 1-52-Buena Vista 1 0
## B 1-52-Huinkul 1 0
## B 1-52-Kennebec 1 0
## B 1-52-S. Rafalela 1 0
## B 25-50 E-B 72-53 A 0 1
## B 25-50 E-Buena Vista 0 1
## B 25-50 E-Huinkul 0 1
## B 25-50 E-Kennebec 0 1
## B 25-50 E-S. Rafalela 0 1
## B 72-53 A-Buena Vista 0 0
## B 72-53 A-Huinkul 0 0
## B 72-53 A-Kennebec 0 0
## B 72-53 A-S. Rafalela 0 0
## Buena Vista-Huinkul 0 0
## Buena Vista-Kennebec 0 0
## Buena Vista-S. Rafalela 0 0
## Huinkul-Kennebec 0 0
## Huinkul-S. Rafalela 0 0
## Kennebec-S. Rafalela 0 0
## variedadeB 72-53 A variedadeBuena Vista
## B 116-51-B 1-52 0 0
## B 116-51-B 25-50 E 0 0
## B 116-51-B 72-53 A -1 0
## B 116-51-Buena Vista 0 -1
## B 116-51-Huinkul 0 0
## B 116-51-Kennebec 0 0
## B 116-51-S. Rafalela 0 0
## B 1-52-B 25-50 E 0 0
## B 1-52-B 72-53 A -1 0
## B 1-52-Buena Vista 0 -1
## B 1-52-Huinkul 0 0
## B 1-52-Kennebec 0 0
## B 1-52-S. Rafalela 0 0
## B 25-50 E-B 72-53 A -1 0
## B 25-50 E-Buena Vista 0 -1
## B 25-50 E-Huinkul 0 0
## B 25-50 E-Kennebec 0 0
## B 25-50 E-S. Rafalela 0 0
## B 72-53 A-Buena Vista 1 -1
## B 72-53 A-Huinkul 1 0
## B 72-53 A-Kennebec 1 0
## B 72-53 A-S. Rafalela 1 0
## Buena Vista-Huinkul 0 1
## Buena Vista-Kennebec 0 1
## Buena Vista-S. Rafalela 0 1
## Huinkul-Kennebec 0 0
## Huinkul-S. Rafalela 0 0
## Kennebec-S. Rafalela 0 0
## variedadeHuinkul variedadeKennebec
## B 116-51-B 1-52 0 0
## B 116-51-B 25-50 E 0 0
## B 116-51-B 72-53 A 0 0
## B 116-51-Buena Vista 0 0
## B 116-51-Huinkul -1 0
## B 116-51-Kennebec 0 -1
## B 116-51-S. Rafalela 0 0
## B 1-52-B 25-50 E 0 0
## B 1-52-B 72-53 A 0 0
## B 1-52-Buena Vista 0 0
## B 1-52-Huinkul -1 0
## B 1-52-Kennebec 0 -1
## B 1-52-S. Rafalela 0 0
## B 25-50 E-B 72-53 A 0 0
## B 25-50 E-Buena Vista 0 0
## B 25-50 E-Huinkul -1 0
## B 25-50 E-Kennebec 0 -1
## B 25-50 E-S. Rafalela 0 0
## B 72-53 A-Buena Vista 0 0
## B 72-53 A-Huinkul -1 0
## B 72-53 A-Kennebec 0 -1
## B 72-53 A-S. Rafalela 0 0
## Buena Vista-Huinkul -1 0
## Buena Vista-Kennebec 0 -1
## Buena Vista-S. Rafalela 0 0
## Huinkul-Kennebec 1 -1
## Huinkul-S. Rafalela 1 0
## Kennebec-S. Rafalela 0 1
## variedadeS. Rafalela
## B 116-51-B 1-52 0
## B 116-51-B 25-50 E 0
## B 116-51-B 72-53 A 0
## B 116-51-Buena Vista 0
## B 116-51-Huinkul 0
## B 116-51-Kennebec 0
## B 116-51-S. Rafalela -1
## B 1-52-B 25-50 E 0
## B 1-52-B 72-53 A 0
## B 1-52-Buena Vista 0
## B 1-52-Huinkul 0
## B 1-52-Kennebec 0
## B 1-52-S. Rafalela -1
## B 25-50 E-B 72-53 A 0
## B 25-50 E-Buena Vista 0
## B 25-50 E-Huinkul 0
## B 25-50 E-Kennebec 0
## B 25-50 E-S. Rafalela -1
## B 72-53 A-Buena Vista 0
## B 72-53 A-Huinkul 0
## B 72-53 A-Kennebec 0
## B 72-53 A-S. Rafalela -1
## Buena Vista-Huinkul 0
## Buena Vista-Kennebec 0
## Buena Vista-S. Rafalela -1
## Huinkul-Kennebec 0
## Huinkul-S. Rafalela -1
## Kennebec-S. Rafalela -1# Passando a matriz de contrastes.
summary(glht(m0, linfct = allctr),
test = adjusted(type = "fdr"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = producao ~ bloco + variedade, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## B 116-51-B 1-52 == 0 0.225 2.067 0.109 0.914357
## B 116-51-B 25-50 E == 0 6.000 2.067 2.903 0.017029 *
## B 116-51-B 72-53 A == 0 -0.300 2.067 -0.145 0.914357
## B 116-51-Buena Vista == 0 10.075 2.067 4.874 0.000251 ***
## B 116-51-Huinkul == 0 -2.550 2.067 -1.234 0.293976
## B 116-51-Kennebec == 0 11.800 2.067 5.708 5.36e-05 ***
## B 116-51-S. Rafalela == 0 -2.950 2.067 -1.427 0.247947
## B 1-52-B 25-50 E == 0 5.775 2.067 2.794 0.019042 *
## B 1-52-B 72-53 A == 0 -0.525 2.067 -0.254 0.898222
## B 1-52-Buena Vista == 0 9.850 2.067 4.765 0.000293 ***
## B 1-52-Huinkul == 0 -2.775 2.067 -1.342 0.271297
## B 1-52-Kennebec == 0 11.575 2.067 5.600 5.91e-05 ***
## B 1-52-S. Rafalela == 0 -3.175 2.067 -1.536 0.216969
## B 25-50 E-B 72-53 A == 0 -6.300 2.067 -3.048 0.013172 *
## B 25-50 E-Buena Vista == 0 4.075 2.067 1.971 0.102125
## B 25-50 E-Huinkul == 0 -8.550 2.067 -4.136 0.001095 **
## B 25-50 E-Kennebec == 0 5.800 2.067 2.806 0.019042 *
## B 25-50 E-S. Rafalela == 0 -8.950 2.067 -4.330 0.000752 ***
## B 72-53 A-Buena Vista == 0 10.375 2.067 5.019 0.000201 ***
## B 72-53 A-Huinkul == 0 -2.250 2.067 -1.088 0.351487
## B 72-53 A-Kennebec == 0 12.100 2.067 5.854 4.62e-05 ***
## B 72-53 A-S. Rafalela == 0 -2.650 2.067 -1.282 0.285098
## Buena Vista-Huinkul == 0 -12.625 2.067 -6.108 3.25e-05 ***
## Buena Vista-Kennebec == 0 1.725 2.067 0.834 0.482294
## Buena Vista-S. Rafalela == 0 -13.025 2.067 -6.301 2.81e-05 ***
## Huinkul-Kennebec == 0 14.350 2.067 6.942 1.04e-05 ***
## Huinkul-S. Rafalela == 0 -0.400 2.067 -0.194 0.913685
## Kennebec-S. Rafalela == 0 -14.750 2.067 -7.136 1.04e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)# Em situações de desbalanceamento ou falta de ortogonalidade nos
# efeitos, é melhor passar a matriz de contrastes do que confiar no
# comportamento default da função. Na dúvida, faça das duas maneiras.
#-----------------------------------------------------------------------
# Representando os resultados em um gráfico de segmentos.
pred <- cbind(lsm$grid, lsm$coef)
# Obtém os intervalos de conbertura individual.
ci <- confint(glht(m0, linfct = lsm$K),
calpha = univariate_calpha())
pred <- cbind(pred, ci$confint)
# Pega as letras.
pred$cld <- cld(g0)$mcletters$Letters
# Faz com que fatores ou characteres tenham os mesmos níveis dos dados
# que alimentaram o modelo.
pred <- equallevels(pred, da)
str(pred)## 'data.frame': 8 obs. of 10 variables:
## $ variedade: Factor w/ 8 levels "B 116-51","B 1-52",..: 1 2 3 4 5 6 7 8
## $ estimate : num 22.5 22.3 16.5 22.8 12.4 ...
## $ se : num 1.46 1.46 1.46 1.46 1.46 ...
## $ df : num 21 21 21 21 21 21 21 21
## $ t.stat : num 15.4 15.2 11.3 15.6 8.5 ...
## $ p.value : num 6.52e-13 7.93e-13 2.23e-10 5.05e-13 3.07e-08 ...
## $ Estimate : num 22.5 22.3 16.5 22.8 12.4 ...
## $ lwr : num 19.46 19.24 13.46 19.76 9.39 ...
## $ upr : num 25.5 25.3 19.5 25.8 15.5 ...
## $ cld : chr "c" "c" "b" "c" ...segplot(reorder(variedade, Estimate) ~ lwr + upr,
centers = Estimate,
data = pred,
xlab = "Produção",
ylab = "Variedade",
draw = FALSE,
cld = pred$cld) +
latticeExtra::layer(
panel.text(x = centers,
y = z,
labels = sprintf("%0.1f %s", centers, cld),
pos = 3))
1.2 Modelo de efeitos aleatórios
library(lme4)
mm0 <- lmer(producao ~ (1 | bloco) + variedade,
data = da)
# Log-verossimilhança.
logLik(m0)## 'log Lik.' -72.99394 (df=12)logLik(mm0)## 'log Lik.' -66.36294 (df=10)# Detalhes do modelo ajustado.
summary(mm0)## Linear mixed model fit by REML ['lmerMod']
## Formula: producao ~ (1 | bloco) + variedade
## Data: da
##
## REML criterion at convergence: 132.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1335 -0.5559 0.1202 0.4370 1.9457
##
## Random effects:
## Groups Name Variance Std.Dev.
## bloco (Intercept) 1.037 1.018
## Residual 8.546 2.923
## Number of obs: 32, groups: bloco, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 22.500 1.548 14.536
## variedadeB 1-52 -0.225 2.067 -0.109
## variedadeB 25-50 E -6.000 2.067 -2.903
## variedadeB 72-53 A 0.300 2.067 0.145
## variedadeBuena Vista -10.075 2.067 -4.874
## variedadeHuinkul 2.550 2.067 1.234
## variedadeKennebec -11.800 2.067 -5.708
## variedadeS. Rafalela 2.950 2.067 1.427
##
## Correlation of Fixed Effects:
## (Intr) vB1-52 vB25-E vB72-A vrddBV vrddHn vrddKn
## variddB1-52 -0.668
## vrddB25-50E -0.668 0.500
## vrddB72-53A -0.668 0.500 0.500
## variddBnVst -0.668 0.500 0.500 0.500
## varieddHnkl -0.668 0.500 0.500 0.500 0.500
## variddKnnbc -0.668 0.500 0.500 0.500 0.500 0.500
## varddS.Rfll -0.668 0.500 0.500 0.500 0.500 0.500 0.500#-----------------------------------------------------------------------
# Médias ajustadas.
# LSmeans(m0, effect = "variedade")
lsm <- LSmeans(mm0, effect = "variedade")
lsm## estimate se df t.stat p.value variedade
## 1 22.500 1.547831 22.18125 14.536469 8.055779e-13 B 116-51
## 2 22.275 1.547831 22.18125 14.391105 9.863334e-13 B 1-52
## 3 16.500 1.547831 22.18125 10.660078 3.398654e-10 B 25-50 E
## 4 22.800 1.547831 22.18125 14.730289 6.165937e-13 B 72-53 A
## 5 12.425 1.547831 22.18125 8.027361 5.256992e-08 Buena Vista
## 6 25.050 1.547831 22.18125 16.183936 9.067232e-14 Huinkul
## 7 10.700 1.547831 22.18125 6.912899 5.827585e-07 Kennebec
## 8 25.450 1.547831 22.18125 16.442362 6.548265e-14 S. Rafalelat(lsm$K)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## (Intercept) 1 1 1 1 1 1 1 1
## variedadeB 1-52 0 1 0 0 0 0 0 0
## variedadeB 25-50 E 0 0 1 0 0 0 0 0
## variedadeB 72-53 A 0 0 0 1 0 0 0 0
## variedadeBuena Vista 0 0 0 0 1 0 0 0
## variedadeHuinkul 0 0 0 0 0 1 0 0
## variedadeKennebec 0 0 0 0 0 0 1 0
## variedadeS. Rafalela 0 0 0 0 0 0 0 1#-----------------------------------------------------------------------
# Contrastes entre médias.
gg0 <- summary(glht(mm0, linfct = mcp(variedade = "Tukey")),
test = adjusted(type = "fdr"))
# Resumo compacto do conjunto de contrastes.
cbind(cld(g0)$mcletters$Letters, cld(gg0)$mcletters$Letters)## [,1] [,2]
## B 116-51 "c" "c"
## B 1-52 "c" "c"
## B 25-50 E "b" "b"
## B 72-53 A "c" "c"
## Buena Vista "ab" "ab"
## Huinkul "c" "c"
## Kennebec "a" "a"
## S. Rafalela "c" "c"#-----------------------------------------------------------------------
# Médias ajustadas.
lsm$K %*% fixef(mm0)## [,1]
## [1,] 22.500
## [2,] 22.275
## [3,] 16.500
## [4,] 22.800
## [5,] 12.425
## [6,] 25.050
## [7,] 10.700
## [8,] 25.450# Matriz de covariância das médias. Covariância surge do amarramento
# dado pelos efeitos aleatórios.
round(lsm$K %*% as.matrix(vcov(mm0)) %*% t(lsm$K), digits = 4)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 2.3958 0.2593 0.2593 0.2593 0.2593 0.2593 0.2593 0.2593
## [2,] 0.2593 2.3958 0.2593 0.2593 0.2593 0.2593 0.2593 0.2593
## [3,] 0.2593 0.2593 2.3958 0.2593 0.2593 0.2593 0.2593 0.2593
## [4,] 0.2593 0.2593 0.2593 2.3958 0.2593 0.2593 0.2593 0.2593
## [5,] 0.2593 0.2593 0.2593 0.2593 2.3958 0.2593 0.2593 0.2593
## [6,] 0.2593 0.2593 0.2593 0.2593 0.2593 2.3958 0.2593 0.2593
## [7,] 0.2593 0.2593 0.2593 0.2593 0.2593 0.2593 2.3958 0.2593
## [8,] 0.2593 0.2593 0.2593 0.2593 0.2593 0.2593 0.2593 2.39582 Delineamento de blocos incompletos balanceados
str(PimentelPg185)## 'data.frame': 30 obs. of 3 variables:
## $ bloc: Factor w/ 10 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ trat: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 1 2 4 1 2 5 1 ...
## $ y : int 35 28 27 30 20 22 28 16 18 36 ...# Diagrama de dispersão.
xyplot(y ~ trat,
groups = bloc, type = "o",
data = PimentelPg185,
auto.key = list(title = "Blocos",
cex.title = 1.1, columns = 5),
xlab = "Tratamento",
ylab = "Resposta")
# Tabela de frequência.
addmargins(xtabs(~trat + bloc, data = PimentelPg185))## bloc
## trat 1 2 3 4 5 6 7 8 9 10 Sum
## 1 1 1 1 1 1 1 0 0 0 0 6
## 2 1 1 1 0 0 0 1 1 1 0 6
## 3 1 0 0 1 1 0 1 1 0 1 6
## 4 0 1 0 1 0 1 1 0 1 1 6
## 5 0 0 1 0 1 1 0 1 1 1 6
## Sum 3 3 3 3 3 3 3 3 3 3 30# Diagrama das conexões entre tratamentos em cada bloco.
k <- nlevels(PimentelPg185$trat)
a <- seq(0, 2 * pi, length.out = k + 1)[-(k + 1)]
par(mfrow = c(2, 5))
col <- 1
for (b in levels(PimentelPg185$bloc)) {
plot(sin(a), cos(a), asp = 1,
xlim = c(-1.1, 1.1),
ylim = c(-1.1, 1.1),
axes = FALSE, xlab = NA, ylab = NA)
mtext(paste("Bloco", b))
i <- unique(as.integer(subset(PimentelPg185, bloc == b)$trat))
cb <- combn(x = i, m = 2)
segments(x0 = sin(a[cb[1, ]]), y0 = cos(a[cb[1, ]]),
x1 = sin(a[cb[2, ]]), y1 = cos(a[cb[2, ]]),
col = col)
text(x = 1.08 * sin(a[i]), y = 1.08 * cos(a[i]),
labels = levels(PimentelPg185$trat)[i])
col <- col + 1
}
layout(1)2.1 Modelo de efeitos fixos
# Tabela de frequência.
pim <- PimentelPg185
str(pim)## 'data.frame': 30 obs. of 3 variables:
## $ bloc: Factor w/ 10 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ trat: Factor w/ 5 levels "1","2","3","4",..: 1 2 3 1 2 4 1 2 5 1 ...
## $ y : int 35 28 27 30 20 22 28 16 18 36 ...m0 <- lm(y ~ bloc + trat, data = pim)
par(mfrow = c(2, 2))
plot(m0)
layout(1)
anova(m0)## Analysis of Variance Table
##
## Response: y
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 9 412.70 45.856 21.383 2.746e-07 ***
## trat 4 283.69 70.922 33.072 1.495e-07 ***
## Residuals 16 34.31 2.144
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)##
## Call:
## lm(formula = y ~ bloc + trat, data = pim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0222 -0.9556 -0.1333 1.1167 1.6889
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.7556 1.0004 35.742 < 2e-16 ***
## bloc2 -5.9111 1.2349 -4.787 0.000202 ***
## bloc3 -9.4444 1.2349 -7.648 9.89e-07 ***
## bloc4 1.3778 1.2349 1.116 0.281022
## bloc5 -7.1556 1.2349 -5.795 2.74e-05 ***
## bloc6 -10.4000 1.2729 -8.170 4.21e-07 ***
## bloc7 0.7778 1.2349 0.630 0.537691
## bloc8 0.2444 1.2349 0.198 0.845578
## bloc9 1.0000 1.2729 0.786 0.443573
## bloc10 -0.3778 1.2729 -0.297 0.770447
## trat2 -9.2000 0.9262 -9.933 3.01e-08 ***
## trat3 -8.0667 0.9262 -8.710 1.81e-07 ***
## trat4 -8.3333 0.9262 -8.998 1.17e-07 ***
## trat5 -7.7333 0.9262 -8.350 3.17e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.464 on 16 degrees of freedom
## Multiple R-squared: 0.953, Adjusted R-squared: 0.9149
## F-statistic: 24.98 on 13 and 16 DF, p-value: 3.72e-08# Médias ajustadas.
lsm <- LSmeans(m0, effect = "trat")
lsm## estimate se df t.stat p.value trat
## 1 32.76667 0.6438886 16 50.88872 3.980184e-19 1
## 2 23.56667 0.6438886 16 36.60053 7.435314e-17 2
## 3 24.70000 0.6438886 16 38.36067 3.535031e-17 3
## 4 24.43333 0.6438886 16 37.94652 4.198038e-17 4
## 5 25.03333 0.6438886 16 38.87836 2.858771e-17 5# Médias amostrais.
aggregate(y ~ trat, data = pim, FUN = mean)## trat y
## 1 1 30.50000
## 2 2 24.33333
## 3 3 26.83333
## 4 4 25.16667
## 5 5 23.66667# Pesos atribuídos aos efeitos de bloco.
t(lsm$K)## [,1] [,2] [,3] [,4] [,5]
## (Intercept) 1.0 1.0 1.0 1.0 1.0
## bloc2 0.1 0.1 0.1 0.1 0.1
## bloc3 0.1 0.1 0.1 0.1 0.1
## bloc4 0.1 0.1 0.1 0.1 0.1
## bloc5 0.1 0.1 0.1 0.1 0.1
## bloc6 0.1 0.1 0.1 0.1 0.1
## bloc7 0.1 0.1 0.1 0.1 0.1
## bloc8 0.1 0.1 0.1 0.1 0.1
## bloc9 0.1 0.1 0.1 0.1 0.1
## bloc10 0.1 0.1 0.1 0.1 0.1
## trat2 0.0 1.0 0.0 0.0 0.0
## trat3 0.0 0.0 1.0 0.0 0.0
## trat4 0.0 0.0 0.0 1.0 0.0
## trat5 0.0 0.0 0.0 0.0 1.0# Médias ajustadas.
lsm$K %*% coef(m0)## [,1]
## [1,] 32.76667
## [2,] 23.56667
## [3,] 24.70000
## [4,] 24.43333
## [5,] 25.03333# Matriz de covariância das médias.
round(lsm$K %*% vcov(m0) %*% t(lsm$K), digits = 4)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.4146 -0.0143 -0.0143 -0.0143 -0.0143
## [2,] -0.0143 0.4146 -0.0143 -0.0143 -0.0143
## [3,] -0.0143 -0.0143 0.4146 -0.0143 -0.0143
## [4,] -0.0143 -0.0143 -0.0143 0.4146 -0.0143
## [5,] -0.0143 -0.0143 -0.0143 -0.0143 0.4146# summary(glht(m0, linfct = mcp(trat = "Tukey")),
# test = adjusted(type = "fdr"))
K <- apc(lsm$K)
summary(glht(m0, linfct = K),
test = adjusted(type = "fdr"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = y ~ bloc + trat, data = pim)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1-2 == 0 9.2000 0.9262 9.933 3.01e-07 ***
## 1-3 == 0 8.0667 0.9262 8.710 6.03e-07 ***
## 1-4 == 0 8.3333 0.9262 8.998 5.85e-07 ***
## 1-5 == 0 7.7333 0.9262 8.350 7.92e-07 ***
## 2-3 == 0 -1.1333 0.9262 -1.224 0.398
## 2-4 == 0 -0.8667 0.9262 -0.936 0.519
## 2-5 == 0 -1.4667 0.9262 -1.584 0.266
## 3-4 == 0 0.2667 0.9262 0.288 0.777
## 3-5 == 0 -0.3333 0.9262 -0.360 0.777
## 4-5 == 0 -0.6000 0.9262 -0.648 0.658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)2.2 Modelo de efeitos aleatórios
mm0 <- lmer(y ~ (1 | bloc) + trat, data = pim)
anova(mm0)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## trat 4 277.86 69.464 32.354summary(mm0)## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ (1 | bloc) + trat
## Data: pim
##
## REML criterion at convergence: 129.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.25326 -0.59067 -0.02004 0.74339 1.10507
##
## Random effects:
## Groups Name Variance Std.Dev.
## bloc (Intercept) 21.114 4.595
## Residual 2.147 1.465
## Number of obs: 30, groups: bloc, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 32.6781 1.5888 20.568
## trat2 -9.0814 0.9237 -9.832
## trat3 -7.8947 0.9237 -8.547
## trat4 -8.2161 0.9237 -8.895
## trat5 -7.6982 0.9237 -8.334
##
## Correlation of Fixed Effects:
## (Intr) trat2 trat3 trat4
## trat2 -0.291
## trat3 -0.291 0.500
## trat4 -0.291 0.500 0.500
## trat5 -0.291 0.500 0.500 0.500# Médias ajustadas.
lsmm <- LSmeans(mm0, effect = "trat")
lsmm## estimate se df t.stat p.value trat
## 1 32.67807 1.589265 11.78748 20.56175 1.345069e-10 1
## 2 23.59663 1.589265 11.78748 14.84751 5.459722e-09 2
## 3 24.78338 1.589265 11.78748 15.59424 3.142449e-09 3
## 4 24.46200 1.589265 11.78748 15.39202 3.640704e-09 4
## 5 24.97992 1.589265 11.78748 15.71791 2.874444e-09 5# Médias ajustadas.
lsmm$K %*% fixef(mm0)## [,1]
## [1,] 32.67807
## [2,] 23.59663
## [3,] 24.78338
## [4,] 24.46200
## [5,] 24.97992# Matriz de covariância das médias. Covariância surge do amarramento
# dado pelos efeitos aleatórios.
round(lsmm$K %*% as.matrix(vcov(mm0)) %*% t(lsmm$K), digits = 4)## [,1] [,2] [,3] [,4] [,5]
## [1,] 2.5242 2.0976 2.0976 2.0976 2.0976
## [2,] 2.0976 2.5242 2.0976 2.0976 2.0976
## [3,] 2.0976 2.0976 2.5242 2.0976 2.0976
## [4,] 2.0976 2.0976 2.0976 2.5242 2.0976
## [5,] 2.0976 2.0976 2.0976 2.0976 2.5242# summary(glht(mm0, linfct = mcp(trat = "Tukey")),
# test = adjusted(type = "fdr"))
K <- apc(lsmm$K)
summary(glht(mm0, linfct = K),
test = adjusted(type = "fdr"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lmer(formula = y ~ (1 | bloc) + trat, data = pim)
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## 1-2 == 0 9.0814 0.9237 9.832 <2e-16 ***
## 1-3 == 0 7.8947 0.9237 8.547 <2e-16 ***
## 1-4 == 0 8.2161 0.9237 8.895 <2e-16 ***
## 1-5 == 0 7.6982 0.9237 8.334 <2e-16 ***
## 2-3 == 0 -1.1867 0.9237 -1.285 0.331
## 2-4 == 0 -0.8654 0.9237 -0.937 0.498
## 2-5 == 0 -1.3833 0.9237 -1.498 0.268
## 3-4 == 0 0.3214 0.9237 0.348 0.809
## 3-5 == 0 -0.1965 0.9237 -0.213 0.832
## 4-5 == 0 -0.5179 0.9237 -0.561 0.719
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)lsm## estimate se df t.stat p.value trat
## 1 32.76667 0.6438886 16 50.88872 3.980184e-19 1
## 2 23.56667 0.6438886 16 36.60053 7.435314e-17 2
## 3 24.70000 0.6438886 16 38.36067 3.535031e-17 3
## 4 24.43333 0.6438886 16 37.94652 4.198038e-17 4
## 5 25.03333 0.6438886 16 38.87836 2.858771e-17 5lsmm## estimate se df t.stat p.value trat
## 1 32.67807 1.589265 11.78748 20.56175 1.345069e-10 1
## 2 23.59663 1.589265 11.78748 14.84751 5.459722e-09 2
## 3 24.78338 1.589265 11.78748 15.59424 3.142449e-09 3
## 4 24.46200 1.589265 11.78748 15.39202 3.640704e-09 4
## 5 24.97992 1.589265 11.78748 15.71791 2.874444e-09 53 Experimento em parcelas subdivididas
- Análise de experimento em parcela subdividida: http://leg.ufpr.br/~walmes/cursoR/geneticaEsalq/script15.html;
- Análise de experimento em parcela subdividida no tempo: http://leg.ufpr.br/~walmes/cursoR/geneticaEsalq/script16.html;
- Análise de VCU (valor de cultivo e uso) de soja: http://leg.ufpr.br/~walmes/cursoR/geneticaEsalq/script19.html;
3.1 Modelo de efeitos fixos
3.2 Modelo de efeitos aleatórios
da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja_planta_tcc.txt",
header = TRUE,
sep = "\t")
str(da)## 'data.frame': 40 obs. of 9 variables:
## $ sistema : Factor w/ 2 levels "convencional",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ adubpk : int 40 40 40 40 60 60 60 60 90 90 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ alt15.12: num 50.6 48.3 51.6 43.3 49 50.3 51 49 48.6 56.3 ...
## $ alt05.01: num 69.6 59.6 68.3 60 63 69.3 64.3 67 58 61.6 ...
## $ alt26.01: num 76.6 74.6 70.3 68 68.3 74 67 64.6 66.3 67.6 ...
## $ alt1vag : num 12.6 14 14 15.3 12.3 13.6 15.3 12.3 15 14 ...
## $ prod : num 2044 2383 3137 2889 2308 ...
## $ p100 : num 32 30.3 33.7 29 29.4 30.9 30.3 31.3 32 30.8 ...xyplot(prod ~ adubpk | bloco,
groups = sistema,
data = da,
type = "o",
auto.key = TRUE)
# Efeito categórico para adub para começar.
da$adub <- factor(da$adubpk)
# Níveis das parcelas.
da$parcela <- with(da, interaction(bloco, sistema, drop = TRUE))
# Adub categórico.
m0 <- lmer(prod ~ bloco + (1 | parcela) + sistema * adub,
data = da)
anova(m0)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## bloco 3 295615 98538 0.6749
## sistema 1 41727 41727 0.2858
## adub 4 1907084 476771 3.2656
## sistema:adub 4 555122 138781 0.9506summary(m0)## Linear mixed model fit by REML ['lmerMod']
## Formula: prod ~ bloco + (1 | parcela) + sistema * adub
## Data: da
##
## REML criterion at convergence: 420.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6868 -0.6002 0.0196 0.6369 2.1150
##
## Random effects:
## Groups Name Variance Std.Dev.
## parcela (Intercept) 60336 245.6
## Residual 145997 382.1
## Number of obs: 40, groups: parcela, 8
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2311.67 291.82 7.922
## blocoII 340.51 299.22 1.138
## blocoIII 371.36 299.22 1.241
## blocoIV 150.04 299.22 0.501
## sistemadireto 85.82 321.20 0.267
## adub60 298.37 270.18 1.104
## adub90 486.35 270.18 1.800
## adub120 354.65 270.18 1.313
## adub150 486.27 270.18 1.800
## sistemadireto:adub60 -206.23 382.10 -0.540
## sistemadireto:adub90 -231.77 382.10 -0.607
## sistemadireto:adub120 387.32 382.10 1.014
## sistemadireto:adub150 187.12 382.10 0.490# Adub linear.
m1 <- update(m0, . ~ bloco + (1 | parcela) + sistema * adubpk)
anova(m0, m1)## Data: da
## Models:
## m1: prod ~ bloco + (1 | parcela) + sistema + adubpk + sistema:adubpk
## m0: prod ~ bloco + (1 | parcela) + sistema * adub
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m1 9 603.04 618.24 -292.52 585.04
## m0 15 611.08 636.41 -290.54 581.08 3.9626 6 0.6817r <- residuals(m0, type = "pearson")
f <- fitted(m0)
b <- unlist(ranef(m0))
grid.arrange(nrow = 2,
xyplot(r ~ f, type = c("p", "smooth")),
xyplot(sqrt(abs(r)) ~ f, type = c("p", "smooth")),
qqmath(r),
qqmath(b))
anova(m0)## Analysis of Variance Table
## Df Sum Sq Mean Sq F value
## bloco 3 295615 98538 0.6749
## sistema 1 41727 41727 0.2858
## adub 4 1907084 476771 3.2656
## sistema:adub 4 555122 138781 0.9506m2 <- update(m0, . ~ bloco + (1 | parcela) + adubpk)
anova(m2, m0)## Data: da
## Models:
## m2: prod ~ bloco + (1 | parcela) + adubpk
## m0: prod ~ bloco + (1 | parcela) + sistema * adub
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## m2 7 601.42 613.24 -293.71 587.42
## m0 15 611.08 636.41 -290.54 581.08 6.3413 8 0.6091pred <- with(da, expand.grid(bloco = levels(bloco),
sistema = levels(sistema),
adubpk = seq(40, 150, by = 10),
KEEP.OUT.ATTRS = FALSE))
str(pred)## 'data.frame': 96 obs. of 3 variables:
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ sistema: Factor w/ 2 levels "convencional",..: 1 1 1 1 2 2 2 2 1 1 ...
## $ adubpk : num 40 40 40 40 40 40 40 40 50 50 ...X <- model.matrix(nobars(formula(m1))[-2], data = pred)
head(X)## (Intercept) blocoII blocoIII blocoIV sistemadireto adubpk
## 1 1 0 0 0 0 40
## 2 1 1 0 0 0 40
## 3 1 0 1 0 0 40
## 4 1 0 0 1 0 40
## 5 1 0 0 0 1 40
## 6 1 1 0 0 1 40
## sistemadireto:adubpk
## 1 0
## 2 0
## 3 0
## 4 0
## 5 40
## 6 40# Para obter o efeito médio dos blocos.
g <- aggregate(X ~ sistema + adubpk, data = pred, FUN = mean)
head(g)## sistema adubpk (Intercept) blocoII blocoIII blocoIV sistemadireto
## 1 convencional 40 1 0.25 0.25 0.25 0
## 2 direto 40 1 0.25 0.25 0.25 1
## 3 convencional 50 1 0.25 0.25 0.25 0
## 4 direto 50 1 0.25 0.25 0.25 1
## 5 convencional 60 1 0.25 0.25 0.25 0
## 6 direto 60 1 0.25 0.25 0.25 1
## adubpk sistemadireto:adubpk
## 1 40 0
## 2 40 40
## 3 50 0
## 4 50 50
## 5 60 0
## 6 60 60w <- 1:2
pred <- g[, w]
X <- as.matrix(g[, -w])
# Verifica se os nomes correspondem.
colnames(X) == names(fixef(m1))## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUEU <- chol(vcov(m1))
pred$se <- sqrt(apply(X %*% t(U), 1, function(x) sum(x^2)))
zval <- qnorm(p = c(lwr = 0.025, fit = 0.5, upr = 0.975))
me <- outer(pred$se, zval, "*")
fit <- crossprod(t(X), fixef(m1))
pred <- cbind(pred, sweep(me, 1, fit, "+"))
str(pred)## 'data.frame': 24 obs. of 6 variables:
## $ sistema: Factor w/ 2 levels "convencional",..: 1 2 1 2 1 2 1 2 1 2 ...
## $ adubpk : num 40 40 50 50 60 60 70 70 80 80 ...
## $ se : num 184 184 173 173 163 ...
## $ lwr : num 2310 2233 2367 2327 2420 ...
## $ fit : num 2670 2593 2705 2665 2740 ...
## $ upr : num 3030 2953 3043 3003 3060 ...sub <-
"Retas não diferem entre si pelo teste de razão de verossimilhanças (5%)."
xyplot(fit ~ adubpk,
groups = sistema,
data = pred,
type = "l",
auto.key = list(points = FALSE,
lines = TRUE,
title = "Sistema de cultivo",
cex.title = 1.1,
corner = c(0.05, 0.95)),
xlab = expression("Adubação com PK" ~ (kg ~ ha^{-1})),
ylab = expression("Produtividade" ~ (kg ~ ha^{-1})),
sub = list(sub, font = 3, cex = 0.8),
prepanel = prepanel.cbH,
cty = "bands",
ly = pred$lwr,
uy = pred$upr,
alpha = 0.1,
fill = 1,
panel.groups = panel.cbH,
panel = panel.superpose)
# labestDataView()
# help(PimentelEg5.2, help_type = "html")
str(RamalhoTb11.1)## 'data.frame': 80 obs. of 3 variables:
## $ bloc: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ...
## $ cult: Factor w/ 16 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ prod: num 10.2 10.7 10.8 12.7 9.3 6.4 10.5 10.6 9.2 5.2 ...addmargins(xtabs(~cult + bloc, data = RamalhoTb11.1))## bloc
## cult 1 2 3 4 Sum
## 1 5 0 0 0 5
## 2 1 4 0 0 5
## 3 1 0 4 0 5
## 4 1 0 0 4 5
## 5 1 2 1 1 5
## 6 1 2 1 1 5
## 7 1 2 1 1 5
## 8 1 2 1 1 5
## 9 1 1 2 1 5
## 10 1 1 2 1 5
## 11 1 1 2 1 5
## 12 1 1 2 1 5
## 13 1 1 1 2 5
## 14 1 1 1 2 5
## 15 1 1 1 2 5
## 16 1 1 1 2 5
## Sum 20 20 20 20 80xyplot(prod ~ cult, data = RamalhoTb11.1,
xlab = "Cultivares de sorgo",
ylab = expression("Produção de grãos"~(kg~parcela^{-1})))
k <- nlevels(RamalhoTb11.1$cult)
a <- seq(0, 2 * pi, length.out = k + 1)[-(k + 1)]
par(mfrow = c(2, 2))
col <- 1
for (b in levels(RamalhoTb11.1$bloc)) {
plot(sin(a), cos(a), asp = 1,
xlim = c(-1.1, 1.1),
ylim = c(-1.1, 1.1),
axes = FALSE, xlab = NA, ylab = NA)
mtext(paste("Bloco", b))
i <- unique(as.integer(subset(RamalhoTb11.1, bloc == b)$cult))
cb <- combn(x = i, m = 2)
segments(x0 = sin(a[cb[1, ]]), y0 = cos(a[cb[1, ]]),
x1 = sin(a[cb[2, ]]), y1 = cos(a[cb[2, ]]),
col = col)
text(x = 1.08 * sin(a[i]), y = 1.08 * cos(a[i]),
labels = levels(RamalhoTb11.1$cult)[i])
col <- col + 1
}
| Extensões de Modelos de Regressão |
Prof. Paulo Justiniano R. Jr & Prof. Walmes M. Zeviani |
