Análise de um experimento em DIC
##-----------------------------------------------------------------------------
## Pacotes.
pkg <- c("lattice", "latticeExtra", "plyr", "reshape",
"doBy", "multcomp")
sapply(pkg, require, character.only=TRUE)
##-----------------------------------------------------------------------------
## Dados de experimento em DIC. Fator ração aplicada à suínos com 4
## níveis. A resposta é o ganho de peso ao final do experimento.
da <-
read.table("http://www.leg.ufpr.br/~walmes/data/pimentel_racoes.txt",
header=TRUE, sep="\t")
str(da)
## 'data.frame': 20 obs. of 2 variables:
## $ racoes : Factor w/ 4 levels "A","B","C","D": 1 1 1 1 1 2 2 2 2 2 ...
## $ ganhopeso: int 35 19 31 15 30 40 35 46 41 33 ...
xyplot(ganhopeso~racoes, data=da, jitter.x=TRUE)
## dotplot(ganhopeso~racoes, data=da, jitter.x=TRUE)
## bwplot(ganhopeso~racoes, data=da, pch="|")
##-----------------------------------------------------------------------------
X <- model.matrix(~racoes, data=da)
## X <- model.matrix(~0+racoes, data=da)
## X <- model.matrix(~-1+racoes, data=da)
## X <- model.matrix(~racoes, data=da, contrasts=list(racoes=contr.SAS))
## X <- model.matrix(~racoes, data=da, contrasts=list(racoes=contr.sum))
## X <- model.matrix(~racoes, data=da, contrasts=list(racoes=contr.helmert))
X
## (Intercept) racoesB racoesC racoesD
## 1 1 0 0 0
## 2 1 0 0 0
## 3 1 0 0 0
## 4 1 0 0 0
## 5 1 0 0 0
## 6 1 1 0 0
## 7 1 1 0 0
## 8 1 1 0 0
## 9 1 1 0 0
## 10 1 1 0 0
## 11 1 0 1 0
## 12 1 0 1 0
## 13 1 0 1 0
## 14 1 0 1 0
## 15 1 0 1 0
## 16 1 0 0 1
## 17 1 0 0 1
## 18 1 0 0 1
## 19 1 0 0 1
## 20 1 0 0 1
## attr(,"assign")
## [1] 0 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$racoes
## [1] "contr.treatment"
## beta = (X'X)^{-1} X'y
XlX <- crossprod(X)
XlX
## (Intercept) racoesB racoesC racoesD
## (Intercept) 20 5 5 5
## racoesB 5 5 0 0
## racoesC 5 0 5 0
## racoesD 5 0 0 5
Xly <- crossprod(X, da$ganhopeso)
Xly
## [,1]
## (Intercept) 595
## racoesB 195
## racoesC 160
## racoesD 110
beta <- solve(XlX, Xly)
beta
## [,1]
## (Intercept) 26
## racoesB 13
## racoesC 6
## racoesD -4
##-----------------------------------------------------------------------------
## Matriz de contrastes aplicados e as estimativas de médias, \mu_i =
## \mu+\tau_i.
Z <- cbind(1, contrasts(C(da$racoes, contr="contr.treatment"))); Z
## B C D
## A 1 0 0 0
## B 1 1 0 0
## C 1 0 1 0
## D 1 0 0 1
Z%*%beta
## [,1]
## A 26
## B 39
## C 32
## D 22
##-----------------------------------------------------------------------------
## Obtendo o quadro de anova matricialmente.
X <- model.matrix(~racoes, data=da)
## Modelo: \mu+\alpha_i
Pmu <- X[,1]%*%solve(crossprod(X[,1]))%*%t(X[,1])
Pmualpha <- X%*%solve(crossprod(X))%*%t(X)
In <- diag(length(da$racoes))
y <- da$ganhopeso
## Somas de quadrados de ração e resídual.
t(y)%*%(Pmualpha-Pmu)%*%y
## [,1]
## [1,] 823.8
t(y)%*%(In-Pmualpha)%*%y
## [,1]
## [1,] 1100
## Graus de liberdade de ração e resídual.
sum(diag(Pmualpha-Pmu))
## [1] 3
sum(diag(In-Pmualpha))
## [1] 16
##-----------------------------------------------------------------------------
## Usando a lm().
m0 <- lm(ganhopeso~racoes, data=da)
coef(m0)
## (Intercept) racoesB racoesC racoesD
## 26 13 6 -4
anova(m0)
## Analysis of Variance Table
##
## Response: ganhopeso
## Df Sum Sq Mean Sq F value Pr(>F)
## racoes 3 824 274.6 3.99 0.027 *
## Residuals 16 1100 68.8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.matrix(m0)
## (Intercept) racoesB racoesC racoesD
## 1 1 0 0 0
## 2 1 0 0 0
## 3 1 0 0 0
## 4 1 0 0 0
## 5 1 0 0 0
## 6 1 1 0 0
## 7 1 1 0 0
## 8 1 1 0 0
## 9 1 1 0 0
## 10 1 1 0 0
## 11 1 0 1 0
## 12 1 0 1 0
## 13 1 0 1 0
## 14 1 0 1 0
## 15 1 0 1 0
## 16 1 0 0 1
## 17 1 0 0 1
## 18 1 0 0 1
## 19 1 0 0 1
## 20 1 0 0 1
## attr(,"assign")
## [1] 0 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$racoes
## [1] "contr.treatment"
##-----------------------------------------------------------------------------
## As médias ajustadas \mu_i = \mu_\alpha_i
LSmeans(m0, effect="racoes")
## estimate se df t.stat p.value lwr upr racoes
## 1 26 3.708 16 7.012 2.935e-06 18.14 33.86 A
## 2 39 3.708 16 10.518 1.355e-08 31.14 46.86 B
## 3 32 3.708 16 8.630 2.047e-07 24.14 39.86 C
## 4 22 3.708 16 5.933 2.103e-05 14.14 29.86 D
LSmatrix(m0, effect="racoes")
## (Intercept) racoesB racoesC racoesD
## [1,] 1 0 0 0
## [2,] 1 1 0 0
## [3,] 1 0 1 0
## [4,] 1 0 0 1
## Matriz de definição dos contrastes ou restrições.
C <- cbind(1, contrasts(da$racoes))
## Matriz de variância e covariância das estimativas: \sigma^2
## (X'X)^{-1}.
V <- vcov(m0)
## Erro padrão das estimativas dos parâmetros.
sqrt(diag(V))
## (Intercept) racoesB racoesC racoesD
## 3.708 5.244 5.244 5.244
## Erro padrão das estimativas das médias.
sqrt(diag(C%*%V%*%t(C)))
## A B C D
## 3.708 3.708 3.708 3.708
##-----------------------------------------------------------------------------
## Estimativa de contrastes.
## \mu_2 - \mu_3.
K <- matrix(c(0,1,-1,0), ncol=1)
t(K)%*%coef(m0)
## [,1]
## [1,] 7
## Erro padrão do contraste.
sqrt(t(K)%*%V%*%K)
## [,1]
## [1,] 5.244
##-----------------------------------------------------------------------------
## Usando a glht() para avaliar hipóteses lineares.
summary(glht(m0, linfct=t(K)))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = ganhopeso ~ racoes, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 7.00 5.24 1.33 0.2
## (Adjusted p values reported -- single-step method)
confint(summary(glht(m0, linfct=C)))
##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = ganhopeso ~ racoes, data = da)
##
## Quantile = 2.783
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## A == 0 26.000 15.682 36.318
## B == 0 39.000 28.682 49.318
## C == 0 32.000 21.682 42.318
## D == 0 22.000 11.682 32.318
## Teste de um conjunto de hipóteses simples.
## \mu_1 - \mu_2, \mu_2 - \mu_3, \mu_3 - \mu_4.
K <- matrix(c(0,-1,0,0,
0,1,-1,0,
0,0,1,-1), ncol=3)
## Aplica uma correção aos p-valores devido a multiplicidade.
summary(glht(m0, linfct=t(K)))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = ganhopeso ~ racoes, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 -13.00 5.24 -2.48 0.064 .
## 2 == 0 7.00 5.24 1.33 0.434
## 3 == 0 10.00 5.24 1.91 0.182
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
## Sem correção, ou seja, p-valor vindo da distribuição t.
summary(glht(m0, linfct=t(K)), test=adjusted(type="none"))
##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = ganhopeso ~ racoes, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 -13.00 5.24 -2.48 0.025 *
## 2 == 0 7.00 5.24 1.33 0.201
## 3 == 0 10.00 5.24 1.91 0.075 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)
## Teste simultâneo de hipótese.
summary(glht(m0, linfct=t(K)), test=Ftest())
##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 1 == 0 -13
## 2 == 0 7
## 3 == 0 10
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 3.99 3 16 0.0267
## Mesmo valor do teste F da anova. Sabe dizer por que?
Análise de um experimento em DBC
##-----------------------------------------------------------------------------
## Experimento feito com a cultura da batatinha, avaliando
## cultivares/variedades. A resposta é a produtividade. Foi instalado em
## DBC.
da <-
read.table("http://www.leg.ufpr.br/~walmes/data/pimentel_batatinha.txt",
header=TRUE, sep="\t")
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",..: 7 6 8 5 3 2 1 4 7 6 ...
## $ producao : num 9.2 21.1 22.6 15.4 12.7 20 23.1 18 13.4 27 ...
## xyplot(producao~variedade, data=da, jitter.x=TRUE)
xyplot(producao~variedade, data=da, jitter.x=TRUE, groups=bloco)
##-----------------------------------------------------------------------------
## Matriz do modelo.
X <- model.matrix(~bloco+variedade, data=da)
head(X)
## (Intercept) blocoII blocoIII blocoIV variedadeB 1-52 variedadeB 25-50 E
## 1 1 0 0 0 0 0
## 2 1 0 0 0 0 0
## 3 1 0 0 0 0 0
## 4 1 0 0 0 0 0
## 5 1 0 0 0 0 1
## 6 1 0 0 0 1 0
## variedadeB 72-53 A variedadeBuena Vista variedadeHuinkul variedadeKennebec
## 1 0 0 0 1
## 2 0 0 1 0
## 3 0 0 0 0
## 4 0 1 0 0
## 5 0 0 0 0
## 6 0 0 0 0
## variedadeS. Rafalela
## 1 0
## 2 0
## 3 1
## 4 0
## 5 0
## 6 0
## beta = (X'X)^{-1} X'y
XlX <- crossprod(X)
XlX
## (Intercept) blocoII blocoIII blocoIV variedadeB 1-52
## (Intercept) 32 8 8 8 4
## blocoII 8 8 0 0 1
## blocoIII 8 0 8 0 1
## blocoIV 8 0 0 8 1
## variedadeB 1-52 4 1 1 1 4
## variedadeB 25-50 E 4 1 1 1 0
## variedadeB 72-53 A 4 1 1 1 0
## variedadeBuena Vista 4 1 1 1 0
## variedadeHuinkul 4 1 1 1 0
## variedadeKennebec 4 1 1 1 0
## variedadeS. Rafalela 4 1 1 1 0
## variedadeB 25-50 E variedadeB 72-53 A variedadeBuena Vista
## (Intercept) 4 4 4
## blocoII 1 1 1
## blocoIII 1 1 1
## blocoIV 1 1 1
## variedadeB 1-52 0 0 0
## variedadeB 25-50 E 4 0 0
## variedadeB 72-53 A 0 4 0
## variedadeBuena Vista 0 0 4
## variedadeHuinkul 0 0 0
## variedadeKennebec 0 0 0
## variedadeS. Rafalela 0 0 0
## variedadeHuinkul variedadeKennebec variedadeS. Rafalela
## (Intercept) 4 4 4
## blocoII 1 1 1
## blocoIII 1 1 1
## blocoIV 1 1 1
## variedadeB 1-52 0 0 0
## variedadeB 25-50 E 0 0 0
## variedadeB 72-53 A 0 0 0
## variedadeBuena Vista 0 0 0
## variedadeHuinkul 4 0 0
## variedadeKennebec 0 4 0
## variedadeS. Rafalela 0 0 4
Xly <- crossprod(X, da$producao)
Xly
## [,1]
## (Intercept) 630.8
## blocoII 170.1
## blocoIII 160.3
## blocoIV 158.3
## variedadeB 1-52 89.1
## variedadeB 25-50 E 66.0
## variedadeB 72-53 A 91.2
## variedadeBuena Vista 49.7
## variedadeHuinkul 100.2
## variedadeKennebec 42.8
## variedadeS. Rafalela 101.8
beta <- solve(XlX, Xly)
##-----------------------------------------------------------------------------
## Quadro de anova.
## Modelo: \mu+\alpha_i+\tau_j.
## \mu: termo sempre presente, modelo nulo (termo 0).
## \alpha's: efeito de bloco (termo 1).
## \tau's: efeito de variedade (termo 2).
a <- attr(X, "assign")
Pmu <- X[,a==0]%*%solve(crossprod(X[,a==0]))%*%t(X[,a==0])
Pmualpha <- X[,a<=1]%*%solve(crossprod(X[,a<=1]))%*%t(X[,a<=1])
Pmualphatau <- X[,a<=2]%*%solve(crossprod(X[,a<=2]))%*%t(X[,a<=2])
In <- diag(length(da$producao))
y <- da$producao
## Somas de quadrados (sequênciais, porém ortogonais por desenho).
t(y)%*%(Pmualphatau-Pmualpha)%*%y ## SQ variedade.
## [,1]
## [1,] 919.7
t(y)%*%(Pmualpha-Pmu)%*%y ## SQ blocos.
## [,1]
## [1,] 50.53
t(y)%*%(In-Pmualphatau)%*%y ## SQ resíduos.
## [,1]
## [1,] 179.5
## É possível fazer uma determinada decompisição matricial de
## Pmualphatau de forma que não seja necessários suas intermediárias.
## Graus de liberdade.
sum(diag(Pmualphatau-Pmualpha)) ## GL variedade.
## [1] 7
sum(diag(Pmualpha-Pmu)) ## GL blocos.
## [1] 3
sum(diag(In-Pmualphatau)) ## GL resíduos.
## [1] 21
##-----------------------------------------------------------------------------
## Usando a função lm().
m0 <- lm(producao~bloco+variedade, data=da)
coef(m0)
## (Intercept) blocoII blocoIII blocoIV
## 20.550 3.500 2.275 2.025
## variedadeB 1-52 variedadeB 25-50 E variedadeB 72-53 A variedadeBuena Vista
## -0.225 -6.000 0.300 -10.075
## variedadeHuinkul variedadeKennebec variedadeS. Rafalela
## 2.550 -11.800 2.950
anova(m0)
## Analysis of Variance Table
##
## Response: producao
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 51 16.8 1.97 0.15
## variedade 7 920 131.4 15.37 5.7e-07 ***
## Residuals 21 179 8.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Médias ajustadas das variedades.
LSmeans(m0, effect="variedade")
## estimate se df t.stat p.value lwr upr variedade
## 1 22.50 1.462 21 15.393 6.524e-13 19.460 25.54 B 116-51
## 2 22.27 1.462 21 15.239 7.925e-13 19.235 25.31 B 1-52
## 3 16.50 1.462 21 11.288 2.229e-10 13.460 19.54 B 25-50 E
## 4 22.80 1.462 21 15.599 5.047e-13 19.760 25.84 B 72-53 A
## 5 12.42 1.462 21 8.501 3.071e-08 9.385 15.46 Buena Vista
## 6 25.05 1.462 21 17.138 8.028e-14 22.010 28.09 Huinkul
## 7 10.70 1.462 21 7.320 3.316e-07 7.660 13.74 Kennebec
## 8 25.45 1.462 21 17.412 5.878e-14 22.410 28.49 S. Rafalela
t(LSmatrix(m0, effect="variedade"))
## [,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