Regressão na presença de cariáveis categóricas
##=============================================================================
## Aplicação de modelos de regressão linear e
## não linear em ciências agrárias
##
## 09 à 11 de Dezembro de 2014 - Goiânia/GO
## Embrapa Arroz e Feijão
##
## Prof. Dr. Walmes M. Zeviani
## LEG/DEST/UFPR
##=============================================================================
##-----------------------------------------------------------------------------
## Definições da sessão.
pkg <- c("lattice", "latticeExtra", "gridExtra", "car", "alr3",
"plyr", "reshape", "doBy", "multcomp", "ellipse", "wzRfun")
sapply(pkg, require, character.only=TRUE)## lattice latticeExtra gridExtra car alr3 plyr
## TRUE TRUE TRUE TRUE TRUE TRUE
## reshape doBy multcomp ellipse wzRfun
## TRUE TRUE TRUE TRUE TRUEtrellis.device(color=FALSE)
extendseq <- function(x, f=0.1, length.out=20){
er <- extendrange(x, f=f)
seq(er[1], er[2], length.out=length.out)
}Regressão linear simples com duas categorias
##-----------------------------------------------------------------------------
## Tamanho de tartarugas em função do sexo.
url <- "http://westfall.ba.ttu.edu/isqs5349/Rdata/turtles.txt"
tur <- read.table(url, header=TRUE, sep="\t")
str(tur)## 'data.frame': 48 obs. of 4 variables:
## $ length: int 98 103 103 105 109 123 123 133 133 133 ...
## $ width : int 81 84 86 86 88 92 95 99 102 102 ...
## $ height: int 38 38 42 42 44 50 46 51 51 51 ...
## $ gender: int 1 1 1 1 1 1 1 1 1 1 ...xtabs(~gender, tur)## gender
## 0 1
## 24 24## xyplot(width+height~length, groups=gender, data=tur, type=c("p","smooth"),
## scales=list(y=list(relation="free")))
splom(tur[,1:3], groups=tur$gender, type=c("p","smooth"), auto.key=TRUE)
##-----------------------------------------------------------------------------
## Ajustar reta separada por gender. São só 2 níveis, então codificada
## como numérica (dummy) ou fator implica no mesmo modelo.
## Possíveis modelos:
## fórmula | n. par. | equação
## ~1 | 1 | b0
## ~length | 2 | b0+b1*x
## ~gender+length | 3 | b0+a0*(sexo==1)+b1*x
## ~gender*length | 4 | b0+a0*(sexo==1)+b1*x+a1*(sexo==1)
## Modelo com interação (modelo maior).
m0 <- lm(height~gender*length, data=tur)
par(mfrow=c(2,2)); plot(m0); layout(1)
summary(m0)##
## Call:
## lm(formula = height ~ gender * length, data = tur)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2913 -0.9333 0.0084 1.0328 3.9903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.17217 3.23946 3.140 0.00302 **
## gender -8.13218 3.89842 -2.086 0.04281 *
## length 0.26934 0.02843 9.475 3.43e-12 ***
## gender:length 0.09821 0.03250 3.022 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.606 on 44 degrees of freedom
## Multiple R-squared: 0.9655, Adjusted R-squared: 0.9631
## F-statistic: 410.5 on 3 and 44 DF, p-value: < 2.2e-16## Usando gender como fator de dois níveis.
tur$Gender <- factor(tur$gender)
str(tur)## 'data.frame': 48 obs. of 5 variables:
## $ length: int 98 103 103 105 109 123 123 133 133 133 ...
## $ width : int 81 84 86 86 88 92 95 99 102 102 ...
## $ height: int 38 38 42 42 44 50 46 51 51 51 ...
## $ gender: int 1 1 1 1 1 1 1 1 1 1 ...
## $ Gender: Factor w/ 2 levels "0","1": 2 2 2 2 2 2 2 2 2 2 ...m1 <- lm(height~Gender*length, data=tur)
m1 <- lm(height~Gender*length, data=tur,
contrasts=list(Gender=contr.SAS))
summary(m1)##
## Call:
## lm(formula = height ~ Gender * length, data = tur, contrasts = list(Gender = contr.SAS))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2913 -0.9333 0.0084 1.0328 3.9903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.03999 2.16877 0.941 0.35204
## Gender1 8.13218 3.89842 2.086 0.04281 *
## length 0.36755 0.01576 23.323 < 2e-16 ***
## Gender1:length -0.09821 0.03250 -3.022 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.606 on 44 degrees of freedom
## Multiple R-squared: 0.9655, Adjusted R-squared: 0.9631
## F-statistic: 410.5 on 3 and 44 DF, p-value: < 2.2e-16## contr.treatment
## contr.SAS
## contr.sum
## contr.helmert
## contr.poly
## Diagnóstico
par(mfrow=c(2,2)); plot(m1); layout(1)
## Essa parametrização de nível de referência é útil para testar
## hipóteses. A parametrização de estimativas por categoria é mais
## interessante para interpretação.
m2 <- update(m1, .~0+Gender/length)
summary(m2)##
## Call:
## lm(formula = height ~ Gender + Gender:length - 1, data = tur,
## contrasts = list(Gender = contr.SAS))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2913 -0.9333 0.0084 1.0328 3.9903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## Gender0 10.17217 3.23946 3.140 0.00302 **
## Gender1 2.03999 2.16877 0.941 0.35204
## Gender0:length 0.26934 0.02843 9.475 3.43e-12 ***
## Gender1:length 0.36755 0.01576 23.323 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.606 on 44 degrees of freedom
## Multiple R-squared: 0.9989, Adjusted R-squared: 0.9988
## F-statistic: 1.031e+04 on 4 and 44 DF, p-value: < 2.2e-16## Veja que o R² aumentou pelo fato de o modelo declarado não ter o
## modelo nulo como modelo mínimo por considerar efeitos iguais a zero.
## Mesmo nessa parametrização é possível comparar as inclinações. Basta
## escrever a matriz de funções lineares corretamente.
L <- rbind(c(0,0,-1,1))
summary(glht(m2, linfct=L))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = height ~ Gender + Gender:length - 1, data = tur,
## contrasts = list(Gender = contr.SAS))
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 0.09821 0.03250 3.022 0.00418 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)##-----------------------------------------------------------------------------
## Predição.
pred <- ddply(tur, .(Gender), summarise, length=extendseq(length, l=12))
pred <- cbind(pred, predict(m2, newdata=pred, interval="confidence"))
p1 <- xyplot(height~length, groups=Gender, data=tur,
xlab="Altura", ylab="Comprimento")
p2 <- xyplot(fit~length, groups=Gender, data=pred, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.5,
prepanel=prepanel.cbH, panel=panel.superpose,
panel.groups=panel.cbH)
p1+as.layer(p2, under=TRUE)
Regressão linear simples com mais de duas categorias
##-----------------------------------------------------------------------------
## Total de tempo dormindo em função do peso corporal e nível de
## agressividade de mamíferos.
str(sleep1)## 'data.frame': 62 obs. of 10 variables:
## $ SWS : num NA 6.3 NA NA 2.1 9.1 15.8 5.2 10.9 8.3 ...
## $ PS : num NA 2 NA NA 1.8 0.7 3.9 1 3.6 1.4 ...
## $ TS : num 3.3 8.3 12.5 16.5 3.9 9.8 19.7 6.2 14.5 9.7 ...
## $ BodyWt : num 6654 1 3.38 0.92 2547 ...
## $ BrainWt: num 5712 6.6 44.5 5.7 4603 ...
## $ Life : num 38.6 4.5 14 NA 69 27 19 30.4 28 50 ...
## $ GP : num 645 42 60 25 624 180 35 392 63 230 ...
## $ P : int 3 3 1 5 3 4 1 4 1 1 ...
## $ SE : int 5 1 1 2 5 4 1 5 2 1 ...
## $ D : int 3 3 1 3 4 4 1 4 1 1 ...## help(sleep1, h="html")
## Danger Index.
sleep1$D <- factor(sleep1$D)
sleep1$lbw <- log(sleep1$BodyWt)
xyplot(TS~lbw|D, data=sleep1)
xyplot(TS~lbw|D, data=sleep1, type=c("p","r"))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lm(TS~D*lbw, data=sleep1)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
## Relação média variância.
## Transformar a resposta?
MASS::boxcox(m0)
m1 <- update(m0, log(TS)~.)
par(mfrow=c(2,2)); plot(m1); layout(1)
## Quadro de anova e de estimativas.
anova(m1)## Analysis of Variance Table
##
## Response: log(TS)
## Df Sum Sq Mean Sq F value Pr(>F)
## D 4 7.6814 1.92036 18.3334 3.348e-09 ***
## lbw 1 2.2637 2.26371 21.6113 2.636e-05 ***
## D:lbw 4 0.5504 0.13761 1.3137 0.2783
## Residuals 48 5.0278 0.10475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m1)##
## Call:
## lm(formula = log(TS) ~ D + lbw + D:lbw, data = sleep1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73527 -0.13732 0.00532 0.21433 0.56142
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.557673 0.084026 30.439 <2e-16 ***
## D2 -0.142271 0.123092 -1.156 0.2535
## D3 -0.326209 0.132434 -2.463 0.0174 *
## D4 -0.310486 0.149163 -2.082 0.0427 *
## D5 -0.734857 0.313919 -2.341 0.0234 *
## lbw -0.041075 0.026081 -1.575 0.1219
## D2:lbw -0.009806 0.068207 -0.144 0.8863
## D3:lbw -0.084594 0.040404 -2.094 0.0416 *
## D4:lbw -0.017844 0.039409 -0.453 0.6527
## D5:lbw -0.079605 0.072602 -1.096 0.2783
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3236 on 48 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.6761, Adjusted R-squared: 0.6154
## F-statistic: 11.13 on 9 and 48 DF, p-value: 3.839e-09## summary(glht(m1))
summary(glht(m1), test=adjusted(type="fdr"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = log(TS) ~ D + lbw + D:lbw, data = sleep1)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) == 0 2.557673 0.084026 30.439 <2e-16 ***
## D2 == 0 -0.142271 0.123092 -1.156 0.3479
## D3 == 0 -0.326209 0.132434 -2.463 0.0781 .
## D4 == 0 -0.310486 0.149163 -2.082 0.0855 .
## D5 == 0 -0.734857 0.313919 -2.341 0.0781 .
## lbw == 0 -0.041075 0.026081 -1.575 0.2031
## D2:lbw == 0 -0.009806 0.068207 -0.144 0.8863
## D3:lbw == 0 -0.084594 0.040404 -2.094 0.0855 .
## D4:lbw == 0 -0.017844 0.039409 -0.453 0.7253
## D5:lbw == 0 -0.079605 0.072602 -1.096 0.3479
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- fdr method)##-----------------------------------------------------------------------------
## Estimativas por grupo.
## R² com SQT corrigida para média.
summary(m1)##
## Call:
## lm(formula = log(TS) ~ D + lbw + D:lbw, data = sleep1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73527 -0.13732 0.00532 0.21433 0.56142
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.557673 0.084026 30.439 <2e-16 ***
## D2 -0.142271 0.123092 -1.156 0.2535
## D3 -0.326209 0.132434 -2.463 0.0174 *
## D4 -0.310486 0.149163 -2.082 0.0427 *
## D5 -0.734857 0.313919 -2.341 0.0234 *
## lbw -0.041075 0.026081 -1.575 0.1219
## D2:lbw -0.009806 0.068207 -0.144 0.8863
## D3:lbw -0.084594 0.040404 -2.094 0.0416 *
## D4:lbw -0.017844 0.039409 -0.453 0.6527
## D5:lbw -0.079605 0.072602 -1.096 0.2783
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3236 on 48 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.6761, Adjusted R-squared: 0.6154
## F-statistic: 11.13 on 9 and 48 DF, p-value: 3.839e-09m2 <- update(m1, .~0+D/lbw)
## R² com SQT **não** corrigida para a média.
summary(m2)##
## Call:
## lm(formula = log(TS) ~ D + D:lbw - 1, data = sleep1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.73527 -0.13732 0.00532 0.21433 0.56142
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## D1 2.55767 0.08403 30.439 < 2e-16 ***
## D2 2.41540 0.08995 26.852 < 2e-16 ***
## D3 2.23146 0.10236 21.799 < 2e-16 ***
## D4 2.24719 0.12324 18.234 < 2e-16 ***
## D5 1.82282 0.30246 6.027 2.28e-07 ***
## D1:lbw -0.04107 0.02608 -1.575 0.121853
## D2:lbw -0.05088 0.06302 -0.807 0.423452
## D3:lbw -0.12567 0.03086 -4.072 0.000173 ***
## D4:lbw -0.05892 0.02954 -1.994 0.051817 .
## D5:lbw -0.12068 0.06776 -1.781 0.081224 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3236 on 48 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.9836, Adjusted R-squared: 0.9802
## F-statistic: 287.5 on 10 and 48 DF, p-value: < 2.2e-16deviance(m1)## [1] 5.027836deviance(m2)## [1] 5.027836m3 <- update(m1, .~0+D+lbw)
summary(m3)##
## Call:
## lm(formula = log(TS) ~ D + lbw - 1, data = sleep1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.76671 -0.17367 0.03659 0.22033 0.59205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## D1 2.59984 0.08007 32.469 < 2e-16 ***
## D2 2.40702 0.08775 27.430 < 2e-16 ***
## D3 2.22816 0.10358 21.512 < 2e-16 ***
## D4 2.27415 0.11370 20.002 < 2e-16 ***
## D5 1.62525 0.13947 11.653 4.07e-16 ***
## lbw -0.07229 0.01574 -4.594 2.80e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3275 on 52 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.9818, Adjusted R-squared: 0.9797
## F-statistic: 467 on 6 and 52 DF, p-value: < 2.2e-16Variáveis indicadoras
##-----------------------------------------------------------------------------
## Preço dos carros em função da categoria (simples, sedan ou cross).
url <- "http://www.leg.ufpr.br/~walmes/data/hb20_venda_webmotors_280314.txt"
db <- read.table(url, header=TRUE, sep="\t")
str(db)## 'data.frame': 783 obs. of 6 variables:
## $ carro : Factor w/ 3 levels "hb20","hb20s",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ especificacao: Factor w/ 20 levels "HYUNDAI HB20 1.0 COMFORT 12V FLEX 4P MANUAL",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ cor : Factor w/ 8 levels "azul","branco",..: 2 2 2 8 2 7 7 8 2 2 ...
## $ km : int 0 0 0 0 0 0 0 0 0 0 ...
## $ anomod : Factor w/ 4 levels "2012/2013","2013/2013",..: 3 3 3 3 2 2 3 4 4 4 ...
## $ preco : int 32890 32900 33240 33240 33250 33250 33295 33490 33890 33990 ...db <- subset(db, select=c("carro", "km", "preco"))
db <- transform(db, km=km/1000, preco=preco/1000)
## Variável indicadora de carro zero km.
## db$novo <- ifelse(db$km==0, 1, 0)
db$novo <- as.integer(db$km==0)
xyplot(preco~km|carro, data=db)
xyplot(preco~km|carro, groups=novo, data=db,
type=c("p", "r"))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
m0 <- lm(preco~carro*(km+novo), data=db)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
summary(m0)##
## Call:
## lm(formula = preco ~ carro * (km + novo), data = db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.6217 -4.1564 -0.6292 4.2707 13.9216
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.11022 0.81051 53.189 < 2e-16 ***
## carrohb20s 5.16347 1.99321 2.591 0.009763 **
## carrohb20x 8.51989 2.12692 4.006 6.78e-05 ***
## km -0.16338 0.04657 -3.508 0.000477 ***
## novo -1.09179 0.85393 -1.279 0.201441
## carrohb20s:km 0.13727 0.16622 0.826 0.409163
## carrohb20x:km -0.01700 0.22733 -0.075 0.940394
## carrohb20s:novo 0.43735 2.05390 0.213 0.831432
## carrohb20x:novo 2.13607 2.23495 0.956 0.339493
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.397 on 774 degrees of freedom
## Multiple R-squared: 0.3358, Adjusted R-squared: 0.329
## F-statistic: 48.92 on 8 and 774 DF, p-value: < 2.2e-16anova(m0)## Analysis of Variance Table
##
## Response: preco
## Df Sum Sq Mean Sq F value Pr(>F)
## carro 2 10887.2 5443.6 186.9200 < 2.2e-16 ***
## km 1 375.3 375.3 12.8884 0.0003515 ***
## novo 1 45.0 45.0 1.5445 0.2143218
## carro:km 2 63.8 31.9 1.0958 0.3347842
## carro:novo 2 26.7 13.4 0.4587 0.6323066
## Residuals 774 22540.9 29.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m1 <- update(m0, .~ carro*km)
summary(m1)##
## Call:
## lm(formula = preco ~ carro + km + carro:km, data = db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7632 -4.2316 -0.5568 4.2432 13.8134
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.12663 0.25499 165.208 < 2e-16 ***
## carrohb20s 5.52513 0.47911 11.532 < 2e-16 ***
## carrohb20x 10.45011 0.65283 16.007 < 2e-16 ***
## km -0.12156 0.03312 -3.670 0.000259 ***
## carrohb20s:km 0.14059 0.09970 1.410 0.158892
## carrohb20x:km -0.12830 0.17801 -0.721 0.471302
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.393 on 777 degrees of freedom
## Multiple R-squared: 0.3341, Adjusted R-squared: 0.3298
## F-statistic: 77.97 on 5 and 777 DF, p-value: < 2.2e-16anova(m1)## Analysis of Variance Table
##
## Response: preco
## Df Sum Sq Mean Sq F value Pr(>F)
## carro 2 10887.2 5443.6 187.1577 < 2.2e-16 ***
## km 1 375.3 375.3 12.9048 0.0003484 ***
## carro:km 2 76.9 38.4 1.3219 0.2672326
## Residuals 777 22599.6 29.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## m2 <- update(m1, .~. -carro:km)
m2 <- update(m1, .~carro+km)
summary(m2)##
## Call:
## lm(formula = preco ~ carro + km, data = db)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.7684 -4.1862 -0.5109 4.1713 13.8413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 42.09873 0.25321 166.263 < 2e-16 ***
## carrohb20s 5.72797 0.45598 12.562 < 2e-16 ***
## carrohb20x 10.35724 0.63180 16.393 < 2e-16 ***
## km -0.11049 0.03077 -3.591 0.00035 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.395 on 779 degrees of freedom
## Multiple R-squared: 0.3318, Adjusted R-squared: 0.3293
## F-statistic: 129 on 3 and 779 DF, p-value: < 2.2e-16##-----------------------------------------------------------------------------
## Predição.
pred <- expand.grid(carro=levels(db$carro),
km=seq(0,50,l=10))
aux <- predict(m2, newdata=pred, interval="confidence")
pred <- cbind(pred, aux)
p1 <- xyplot(preco~km, groups=carro, data=db, auto.key=TRUE)
p2 <- xyplot(fit~km, groups=carro, data=pred, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.5,
panel=panel.superpose, panel.groups=panel.cbH)
p1+as.layer(p2, under=TRUE)
Regressão na análise de experimentos
Exemplo 1
##-----------------------------------------------------------------------------
## Dados de índice agronômico de milho.
da <- read.table("http://www.leg.ufpr.br/~walmes/cursoR/anovareg.txt",
header=TRUE, sep="\t")
str(da)## 'data.frame': 72 obs. of 4 variables:
## $ cultivar: Factor w/ 3 levels "Ag-1002","BR-300",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dose : int 0 0 0 0 60 60 60 60 120 120 ...
## $ bloco : Factor w/ 4 levels "I","II","III",..: 1 2 3 4 1 2 3 4 1 2 ...
## $ indice : int 46 48 44 46 48 47 49 48 52 50 ...xyplot(indice~dose|cultivar, data=da, type=c("p","a"), jitter.x=TRUE)
## xyplot(indice~dose|cultivar, data=da, groups=bloco, type="b")
## Por que sempre dois pontos são iguais?
##-----------------------------------------------------------------------------
## Ajuste do modelo saturado (equivalente ao poli grau 5).
da$Dose <- factor(da$dose)
m0 <- lm(indice~bloco+cultivar*Dose, data=da)
## Diagnóstico.
par(mfrow=c(2,2)); plot(m0); layout(1)
anova(m0)## Analysis of Variance Table
##
## Response: indice
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 13.83 4.611 4.4232 0.007709 **
## cultivar 2 38.08 19.042 18.2657 1.042e-06 ***
## Dose 5 729.17 145.833 139.8903 < 2.2e-16 ***
## cultivar:Dose 10 53.25 5.325 5.1080 3.863e-05 ***
## Residuals 51 53.17 1.042
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Ajuste dos modelos de trabalho.
m1 <- update(m0, .~bloco+cultivar*(dose+I(dose^2)+I(dose^3)))
m2 <- update(m0, .~bloco+cultivar*(dose+I(dose^2)))
## Teste da falta de ajuste.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: indice ~ bloco + cultivar * Dose
## Model 2: indice ~ bloco + cultivar + dose + I(dose^2) + I(dose^3) + cultivar:dose +
## cultivar:I(dose^2) + cultivar:I(dose^3)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 51 53.167
## 2 57 57.542 -6 -4.375 0.6995 0.6512anova(m0, m2)## Analysis of Variance Table
##
## Model 1: indice ~ bloco + cultivar * Dose
## Model 2: indice ~ bloco + cultivar + dose + I(dose^2) + cultivar:dose +
## cultivar:I(dose^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 51 53.167
## 2 60 64.804 -9 -11.637 1.2404 0.2922##-----------------------------------------------------------------------------
## Predição (considerar o efeito de um bloco qualquer, o 1 primeiro).
pred <- with(da,
expand.grid(bloco=levels(bloco)[1],
cultivar=levels(cultivar),
dose=extendseq(dose)))
pred <- cbind(pred, predict(m2, newdata=pred, interval="confidence"))
p1 <- xyplot(indice~dose|cultivar, data=da)
p2 <- xyplot(fit~dose|cultivar, data=pred, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.25,
prepanel=prepanel.cbH,
panel=panel.cbH)
p1+as.layer(p2)
## Não passa no meio dos pontos porque essa predição está
## particularizada para o efeito do bloco 1 e não na média dos blocos.
##-----------------------------------------------------------------------------
## Caso alguém tenha necessidade, é possível conduzir testes para saber
## se as inclinações são iguais na origem, ou seja, comparar o b1 das
## curvas. Só é necessário montar a matriz de contrastes.
levels(da$cultivar)## [1] "Ag-1002" "BR-300" "Pioneer-B815"cbind(coef(m2))## [,1]
## (Intercept) 4.658929e+01
## blocoII -1.055556e+00
## blocoIII -8.888889e-01
## blocoIV -1.055556e+00
## cultivarBR-300 -3.491071e+00
## cultivarPioneer-B815 1.339286e+00
## dose 4.336310e-02
## I(dose^2) -5.208333e-05
## cultivarBR-300:dose 5.132440e-02
## cultivarPioneer-B815:dose 1.223214e-02
## cultivarBR-300:I(dose^2) -1.401290e-04
## cultivarPioneer-B815:I(dose^2) -5.704365e-05L <- rbind("Ag vs BR"=c(0, 0,0,0, 0,0, 0, 0, 1,0, 0,0),
"Ag vs Pi"=c(0, 0,0,0, 0,0, 0, 0, 0,1, 0,0),
"Br vs Pi"=c(0, 0,0,0, 0,0, 0, 0,-1,1, 0,0))
## Estimativa das diferenças em inclinação na origem.
L%*%coef(m2)## [,1]
## Ag vs BR 0.05132440
## Ag vs Pi 0.01223214
## Br vs Pi -0.03909226summary(glht(m2, linfct=L))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = indice ~ bloco + cultivar + dose + I(dose^2) + cultivar:dose +
## cultivar:I(dose^2), data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Ag vs BR == 0 0.05132 0.01044 4.915 < 1e-04 ***
## Ag vs Pi == 0 0.01223 0.01044 1.171 0.47467
## Br vs Pi == 0 -0.03909 0.01044 -3.744 0.00114 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)## A produtividade na dose 240 é a mesma nessas cultivares? Como fazer
## esse teste de médias?
LSmeans(m2, effect="cultivar",
at=list("dose"=240))## estimate se df t.stat p.value cultivar dose I(dose^2)
## 1 53.24643 0.2879825 60 184.8947 2.078013e-84 Ag-1002 240 33000
## 2 54.00179 0.2879825 60 187.5176 8.937457e-85 BR-300 240 33000
## 3 54.23571 0.2879825 60 188.3299 6.898701e-85 Pioneer-B815 240 33000L <- LSmatrix(m2, effect="cultivar",
at=list("dose"=240))
levels(da$cultivar)## [1] "Ag-1002" "BR-300" "Pioneer-B815"L <- rbind(
"Ag"=LSmatrix(m2, at=list("cultivar"="Ag-1002", "dose"=320)),
"Br"=LSmatrix(m2, at=list("cultivar"="BR-300", "dose"=250)),
"Pi"=LSmatrix(m2, at=list("cultivar"="Pioneer-B815", "dose"=255)))
## Estimativas das produtividades nessa dose.
summary(glht(m2, linfct=L))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = indice ~ bloco + cultivar + dose + I(dose^2) + cultivar:dose +
## cultivar:I(dose^2), data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 54.3821 0.5982 90.91 <2e-16 ***
## 2 == 0 54.0068 0.2967 182.06 <2e-16 ***
## 3 == 0 54.2594 0.3040 178.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)## Estimativas dos contrastes.
summary(glht(m2, linfct=apc(L)))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = indice ~ bloco + cultivar + dose + I(dose^2) + cultivar:dose +
## cultivar:I(dose^2), data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1-2 == 0 0.3753 0.6677 0.562 0.837
## 1-3 == 0 0.1228 0.6710 0.183 0.981
## 2-3 == 0 -0.2526 0.4247 -0.595 0.820
## (Adjusted p values reported -- single-step method)summary(glht(m2, linfct=apc(L)), test=Ftest())##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 1-2 == 0 0.3753
## 1-3 == 0 0.1228
## 2-3 == 0 -0.2526
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 0.2557 2 60 0.7752##-----------------------------------------------------------------------------
## Como obter as equações (se é que fato são necessárias).
m3 <- update(m2, .~0+cultivar/(dose+I(dose^2))+bloco,
contrasts=list(bloco=contr.sum))
summary(m3)$coeff## Estimate Std. Error t value Pr(>|t|)
## cultivarAg-1002 4.583929e+01 4.709564e-01 97.3323244 9.533156e-68
## cultivarBR-300 4.234821e+01 4.709564e-01 89.9195977 1.070432e-65
## cultivarPioneer-B815 4.717857e+01 4.709564e-01 100.1760815 1.711340e-68
## bloco1 7.500000e-01 2.121389e-01 3.5354202 7.914710e-04
## bloco2 -3.055556e-01 2.121389e-01 -1.4403564 1.549646e-01
## bloco3 -1.388889e-01 2.121389e-01 -0.6547075 5.151582e-01
## cultivarAg-1002:dose 4.336310e-02 7.383254e-03 5.8731689 1.999069e-07
## cultivarBR-300:dose 9.468750e-02 7.383254e-03 12.8246306 7.687617e-19
## cultivarPioneer-B815:dose 5.559524e-02 7.383254e-03 7.5299104 3.112663e-10
## cultivarAg-1002:I(dose^2) -5.208333e-05 2.362354e-05 -2.2047219 3.131888e-02
## cultivarBR-300:I(dose^2) -1.922123e-04 2.362354e-05 -8.1364736 2.867760e-11
## cultivarPioneer-B815:I(dose^2) -1.091270e-04 2.362354e-05 -4.6194173 2.085883e-05split(coef(m3), m3$assign)## $`1`
## cultivarAg-1002 cultivarBR-300 cultivarPioneer-B815
## 45.83929 42.34821 47.17857
##
## $`2`
## bloco1 bloco2 bloco3
## 0.7500000 -0.3055556 -0.1388889
##
## $`3`
## cultivarAg-1002:dose cultivarBR-300:dose cultivarPioneer-B815:dose
## 0.04336310 0.09468750 0.05559524
##
## $`4`
## cultivarAg-1002:I(dose^2) cultivarBR-300:I(dose^2)
## -5.208333e-05 -1.922123e-04
## cultivarPioneer-B815:I(dose^2)
## -1.091270e-04##-----------------------------------------------------------------------------
## É possível também, caso considerar necessário, repartir as somas de
## quadrados.
m4 <- aov(indice~bloco+cultivar/(dose+I(dose^2)), data=da)
m4$assign## [1] 0 1 1 1 2 2 3 3 3 4 4 4anova(m4)## Analysis of Variance Table
##
## Response: indice
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 13.83 4.611 4.2693 0.008472 **
## cultivar 2 38.08 19.042 17.6300 9.489e-07 ***
## cultivar:dose 3 670.98 223.660 207.0788 < 2.2e-16 ***
## cultivar:I(dose^2) 3 99.80 33.267 30.8007 3.524e-12 ***
## Residuals 60 64.80 1.080
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coef(m4)## (Intercept) blocoII
## 4.658929e+01 -1.055556e+00
## blocoIII blocoIV
## -8.888889e-01 -1.055556e+00
## cultivarBR-300 cultivarPioneer-B815
## -3.491071e+00 1.339286e+00
## cultivarAg-1002:dose cultivarBR-300:dose
## 4.336310e-02 9.468750e-02
## cultivarPioneer-B815:dose cultivarAg-1002:I(dose^2)
## 5.559524e-02 -5.208333e-05
## cultivarBR-300:I(dose^2) cultivarPioneer-B815:I(dose^2)
## -1.922123e-04 -1.091270e-04summary(m4,
split=list(
"cultivar:dose"=list(Ag=1, BR=2, Pi=3),
"cultivar:I(dose^2)"=list(BR=2, Pi=3, Ag=1)
))## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 13.8 4.6 4.269 0.00847 **
## cultivar 2 38.1 19.0 17.630 9.49e-07 ***
## cultivar:dose 3 671.0 223.7 207.079 < 2e-16 ***
## cultivar:dose: Ag 1 193.9 193.9 179.516 < 2e-16 ***
## cultivar:dose: BR 1 345.4 345.4 319.824 < 2e-16 ***
## cultivar:dose: Pi 1 131.7 131.7 121.897 4.41e-16 ***
## cultivar:I(dose^2) 3 99.8 33.3 30.801 3.52e-12 ***
## cultivar:I(dose^2): BR 1 71.5 71.5 66.202 2.87e-11 ***
## cultivar:I(dose^2): Pi 1 23.0 23.0 21.339 2.09e-05 ***
## cultivar:I(dose^2): Ag 1 5.2 5.2 4.861 0.03132 *
## Residuals 60 64.8 1.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m5 <- aov(indice~bloco+cultivar/(dose+I(dose^2)+I(dose^3)), data=da)
summary(m5,
split=list(
"cultivar:dose"=list(Ag=1, BR=2, Pi=3),
"cultivar:I(dose^2)"=list(Ag=1, BR=2, Pi=3),
"cultivar:I(dose^3)"=list(Ag=1, BR=2, Pi=3)
))## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 3 13.8 4.6 4.568 0.00617 **
## cultivar 2 38.1 19.0 18.862 5.17e-07 ***
## cultivar:dose 3 671.0 223.7 221.554 < 2e-16 ***
## cultivar:dose: Ag 1 193.9 193.9 192.064 < 2e-16 ***
## cultivar:dose: BR 1 345.4 345.4 342.180 < 2e-16 ***
## cultivar:dose: Pi 1 131.7 131.7 130.418 2.32e-16 ***
## cultivar:I(dose^2) 3 99.8 33.3 32.954 1.74e-12 ***
## cultivar:I(dose^2): Ag 1 5.2 5.2 5.201 0.02634 *
## cultivar:I(dose^2): BR 1 71.5 71.5 70.830 1.41e-11 ***
## cultivar:I(dose^2): Pi 1 23.0 23.0 22.831 1.28e-05 ***
## cultivar:I(dose^3) 3 7.3 2.4 2.398 0.07740 .
## cultivar:I(dose^3): Ag 1 0.7 0.7 0.666 0.41788
## cultivar:I(dose^3): BR 1 0.2 0.2 0.166 0.68479
## cultivar:I(dose^3): Pi 1 6.4 6.4 6.362 0.01448 *
## Residuals 57 57.5 1.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1levels(da$cultivar)## [1] "Ag-1002" "BR-300" "Pioneer-B815"da$ind <- as.integer(da$cultivar=="Pioneer-B815")
m6 <- lm(indice~bloco+cultivar*(dose+I(dose^2))+ind:I(dose^3), data=da)
summary(m6)##
## Call:
## lm(formula = indice ~ bloco + cultivar * (dose + I(dose^2)) +
## ind:I(dose^3), data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.08750 -0.59226 -0.02361 0.34861 2.96806
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.659e+01 4.944e-01 94.233 < 2e-16 ***
## blocoII -1.056e+00 3.316e-01 -3.183 0.00232 **
## blocoIII -8.889e-01 3.316e-01 -2.681 0.00951 **
## blocoIV -1.056e+00 3.316e-01 -3.183 0.00232 **
## cultivarBR-300 -3.491e+00 6.375e-01 -5.476 9.39e-07 ***
## cultivarPioneer-B815 1.812e+00 6.639e-01 2.729 0.00837 **
## dose 4.336e-02 7.067e-03 6.136 7.68e-08 ***
## I(dose^2) -5.208e-05 2.261e-05 -2.303 0.02480 *
## cultivarBR-300:dose 5.132e-02 9.994e-03 5.135 3.32e-06 ***
## cultivarPioneer-B815:dose -2.371e-02 1.729e-02 -1.371 0.17547
## cultivarBR-300:I(dose^2) -1.401e-04 3.198e-05 -4.382 4.90e-05 ***
## cultivarPioneer-B815:I(dose^2) 2.709e-04 1.326e-04 2.042 0.04559 *
## ind:I(dose^3) -7.287e-07 2.861e-07 -2.548 0.01347 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9947 on 59 degrees of freedom
## Multiple R-squared: 0.9342, Adjusted R-squared: 0.9208
## F-statistic: 69.82 on 12 and 59 DF, p-value: < 2.2e-16## As mesmas conclusões são obtidas se olhar para o quadro de
## estimativas.
summary.lm(m4)$coeff## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.658929e+01 5.165296e-01 90.196732 8.911513e-66
## blocoII -1.055556e+00 3.464213e-01 -3.047028 3.433041e-03
## blocoIII -8.888889e-01 3.464213e-01 -2.565919 1.280581e-02
## blocoIV -1.055556e+00 3.464213e-01 -3.047028 3.433041e-03
## cultivarBR-300 -3.491071e+00 6.660330e-01 -5.241589 2.170142e-06
## cultivarPioneer-B815 1.339286e+00 6.660330e-01 2.010840 4.884237e-02
## cultivarAg-1002:dose 4.336310e-02 7.383254e-03 5.873169 1.999069e-07
## cultivarBR-300:dose 9.468750e-02 7.383254e-03 12.824631 7.687617e-19
## cultivarPioneer-B815:dose 5.559524e-02 7.383254e-03 7.529910 3.112663e-10
## cultivarAg-1002:I(dose^2) -5.208333e-05 2.362354e-05 -2.204722 3.131888e-02
## cultivarBR-300:I(dose^2) -1.922123e-04 2.362354e-05 -8.136474 2.867760e-11
## cultivarPioneer-B815:I(dose^2) -1.091270e-04 2.362354e-05 -4.619417 2.085883e-05##-----------------------------------------------------------------------------
## Obter os valores preditos na média dos blocos.
pred <- with(da,
expand.grid(bloco=levels(bloco)[1],
cultivar=levels(cultivar),
dose=extendseq(dose)))
pred$ind <- as.integer(pred$cultivar=="Pioneer-B815")
str(pred)## 'data.frame': 60 obs. of 4 variables:
## $ bloco : Factor w/ 1 level "I": 1 1 1 1 1 1 1 1 1 1 ...
## $ cultivar: Factor w/ 3 levels "Ag-1002","BR-300",..: 1 2 3 1 2 3 1 2 3 1 ...
## $ dose : num -30 -30 -30 -11.1 -11.1 ...
## $ ind : int 0 0 1 0 0 1 0 0 1 0 ...
## - attr(*, "out.attrs")=List of 2
## ..$ dim : Named int 1 3 20
## .. ..- attr(*, "names")= chr "bloco" "cultivar" "dose"
## ..$ dimnames:List of 3
## .. ..$ bloco : chr "bloco=I"
## .. ..$ cultivar: chr "cultivar=Ag-1002" "cultivar=BR-300" "cultivar=Pioneer-B815"
## .. ..$ dose : chr "dose=-30.000000" "dose=-11.052632" "dose= 7.894737" "dose= 26.842105" ...pred <- equallevels(pred, da)
m2 <- m6
L <- model.matrix(formula(m2)[-2], data=pred)
head(L)## (Intercept) blocoII blocoIII blocoIV cultivarBR-300 cultivarPioneer-B815 dose
## 1 1 0 0 0 0 0 -30.00000
## 2 1 0 0 0 1 0 -30.00000
## 3 1 0 0 0 0 1 -30.00000
## 4 1 0 0 0 0 0 -11.05263
## 5 1 0 0 0 1 0 -11.05263
## 6 1 0 0 0 0 1 -11.05263
## I(dose^2) cultivarBR-300:dose cultivarPioneer-B815:dose cultivarBR-300:I(dose^2)
## 1 900.0000 0.00000 0.00000 0.0000
## 2 900.0000 -30.00000 0.00000 900.0000
## 3 900.0000 0.00000 -30.00000 0.0000
## 4 122.1607 0.00000 0.00000 0.0000
## 5 122.1607 -11.05263 0.00000 122.1607
## 6 122.1607 0.00000 -11.05263 0.0000
## cultivarPioneer-B815:I(dose^2) ind:I(dose^3)
## 1 0.0000 0.000
## 2 0.0000 0.000
## 3 900.0000 -27000.000
## 4 0.0000 0.000
## 5 0.0000 0.000
## 6 122.1607 -1350.197## Nas colunas de correspondem ao efeito dos blocos, colocar 0.25 (1/4).
m2$assign ## Coeficientes de bloco são indicados por 1.## [1] 0 1 1 1 2 2 3 4 5 5 6 6 7L[,m2$assign==1] <- L[,m2$assign==1]+0.25
head(L)## (Intercept) blocoII blocoIII blocoIV cultivarBR-300 cultivarPioneer-B815 dose
## 1 1 0.25 0.25 0.25 0 0 -30.00000
## 2 1 0.25 0.25 0.25 1 0 -30.00000
## 3 1 0.25 0.25 0.25 0 1 -30.00000
## 4 1 0.25 0.25 0.25 0 0 -11.05263
## 5 1 0.25 0.25 0.25 1 0 -11.05263
## 6 1 0.25 0.25 0.25 0 1 -11.05263
## I(dose^2) cultivarBR-300:dose cultivarPioneer-B815:dose cultivarBR-300:I(dose^2)
## 1 900.0000 0.00000 0.00000 0.0000
## 2 900.0000 -30.00000 0.00000 900.0000
## 3 900.0000 0.00000 -30.00000 0.0000
## 4 122.1607 0.00000 0.00000 0.0000
## 5 122.1607 -11.05263 0.00000 122.1607
## 6 122.1607 0.00000 -11.05263 0.0000
## cultivarPioneer-B815:I(dose^2) ind:I(dose^3)
## 1 0.0000 0.000
## 2 0.0000 0.000
## 3 900.0000 -27000.000
## 4 0.0000 0.000
## 5 0.0000 0.000
## 6 122.1607 -1350.197## Obter os ic com a glht.
ic <- confint(glht(m2, linfct=L), calpha=univariate_calpha())$confint
pred <- cbind(pred, ic)
p3 <- xyplot(Estimate~dose|cultivar, data=pred, type="l",
ly=pred$lwr, uy=pred$upr, cty="bands",
alpha=0.25, col=2, fill=2,
prepanel=prepanel.cbH,
panel=panel.cbH)
p1+as.layer(p2)+as.layer(p3)
## As bandas de confiança no efeito ponderado dos blocos são visualmente
## apelativas pois são mais precisas e passam entre os pontos.Exemplo 2
##-----------------------------------------------------------------------------
## Dados de experimento com soja sob 3 níveis de umidade do solo (massa
## de água/massa de solo) e níveis de potássio fornecidos na adubação. A
## hipótese sob estudo é que presença de potássio dá suporte para a
## cultura superar a ocorrência de períodos sem chuva. Experimento
## conduzido em casa de vegetação com 5 blocos, 2 plantas por
## u.e. (vaso).
## Para acessar o artigo (pdf online).
## browseURL("http://www.scielo.br/pdf/rca/v43n2/a03v43n2.pdf")
da <- read.table("http://www.leg.ufpr.br/~walmes/data/soja.txt",
header=TRUE, sep="\t", dec=",")
names(da) <- substr(names(da), 1, 4)
str(da)## 'data.frame': 75 obs. of 10 variables:
## $ pota: int 0 30 60 120 180 0 30 60 120 180 ...
## $ agua: num 37.5 37.5 37.5 37.5 37.5 50 50 50 50 50 ...
## $ bloc: Factor w/ 5 levels "I","II","III",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ reng: num 14.6 21.5 24.6 21.9 28.1 ...
## $ peso: num 10.7 13.5 15.8 12.8 14.8 ...
## $ kgra: num 15.1 17.1 19.1 18.1 19.1 ...
## $ pgra: num 1.18 0.99 0.82 0.85 0.88 1.05 1.08 0.74 1.01 1.01 ...
## $ ts : int 136 159 156 171 190 140 193 200 208 237 ...
## $ nvi : int 22 2 0 2 0 20 6 6 7 10 ...
## $ nv : int 56 62 66 68 82 63 86 94 86 97 ...xyplot(reng~pota, groups=agua, data=da,
type=c("p","a"))
xyplot(reng~pota|agua, data=da,
type=c("p","a"))
##-----------------------------------------------------------------------------
## Ajuste do modelo.
da$A <- factor(da$agua)
da$K <- factor(da$pota)
## Ajuste.
m0 <- lm(reng~bloc+A*K, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Remove obs 74 (uma planta que definou).
da <- da[-74,]
m0 <- update(m0, data=da)
par(mfrow=c(2,2)); plot(m0); layout(1)
## Quadro de anova.
anova(m0)## Analysis of Variance Table
##
## Response: reng
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 100.39 25.10 4.8948 0.001896 **
## A 2 446.53 223.26 43.5456 4.617e-12 ***
## K 4 2592.13 648.03 126.3925 < 2.2e-16 ***
## A:K 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## Estimativas dos parâmentros (efeitos).
summary(m0)##
## Call:
## lm(formula = reng ~ bloc + A * K, data = da)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4674 -1.4230 0.3613 1.2828 5.3249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.0727 1.1402 12.342 < 2e-16 ***
## blocII 1.0660 0.8268 1.289 0.202694
## blocIII -0.8607 0.8268 -1.041 0.302455
## blocIV -1.4433 0.8268 -1.746 0.086454 .
## blocV -1.5256 0.8451 -1.805 0.076508 .
## A50 2.1920 1.4321 1.531 0.131591
## A62.5 2.5700 1.4321 1.795 0.078215 .
## K30 6.8140 1.4321 4.758 1.46e-05 ***
## K60 10.4060 1.4321 7.266 1.38e-09 ***
## K120 11.7880 1.4321 8.231 3.67e-11 ***
## K180 11.8700 1.4321 8.289 2.96e-11 ***
## A50:K30 0.0440 2.0253 0.022 0.982746
## A62.5:K30 -1.9160 2.0253 -0.946 0.348263
## A50:K60 2.5740 2.0253 1.271 0.209099
## A62.5:K60 3.3320 2.0253 1.645 0.105630
## A50:K120 2.2860 2.0253 1.129 0.263907
## A62.5:K120 8.4613 2.0920 4.045 0.000165 ***
## A50:K180 1.1780 2.0253 0.582 0.563178
## A62.5:K180 9.1860 2.0253 4.536 3.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.264 on 55 degrees of freedom
## Multiple R-squared: 0.9239, Adjusted R-squared: 0.899
## F-statistic: 37.11 on 18 and 55 DF, p-value: < 2.2e-16##-----------------------------------------------------------------------------
## Polinômio de segundo grau para o efeito de potássio.
## m1 <- lm(reng~bloc+A*(pota+I(pota^2)), data=da)
m1 <- update(m0, .~bloc+A*(pota+I(pota^2)))
par(mfrow=c(2,2)); plot(m1); layout(1)
## Testa o abandono dos termos (falta de ajuste).
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: reng ~ bloc + A * K
## Model 2: reng ~ bloc + A + pota + I(pota^2) + A:pota + A:I(pota^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 55 281.99
## 2 61 352.46 -6 -70.47 2.2907 0.04804 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m1)## Analysis of Variance Table
##
## Response: reng
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 100.39 25.10 4.3434 0.003718 **
## A 2 446.53 223.26 38.6399 1.443e-11 ***
## pota 1 1999.89 1999.89 346.1175 < 2.2e-16 ***
## I(pota^2) 1 555.07 555.07 96.0640 3.821e-14 ***
## A:pota 2 245.17 122.58 21.2152 1.013e-07 ***
## A:I(pota^2) 2 7.53 3.77 0.6516 0.524783
## Residuals 61 352.46 5.78
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Com b2 comum aos níveis de água.
## m2 <- lm(reng~bloc+A*pota+I(pota^2), data=da)
m2 <- update(m1, .~bloc+A*pota+I(pota^2))
anova(m2, m0)## Analysis of Variance Table
##
## Model 1: reng ~ bloc + A + pota + I(pota^2) + A:pota
## Model 2: reng ~ bloc + A * K
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 63 359.99
## 2 55 281.99 8 78 1.9017 0.07823 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Predição com bandas (considerando efeito de bloco I).
pred <- expand.grid(bloc="I",
A=levels(da$A),
pota=seq(0, 180, by=3))
aux <- predict(m2, newdata=pred, interval="confidence")
pred <- cbind(pred, aux)
str(pred)## 'data.frame': 183 obs. of 6 variables:
## $ bloc: Factor w/ 1 level "I": 1 1 1 1 1 1 1 1 1 1 ...
## $ A : Factor w/ 3 levels "37.5","50","62.5": 1 2 3 1 2 3 1 2 3 1 ...
## $ pota: num 0 0 0 3 3 3 6 6 6 9 ...
## $ fit : num 14.4 17.2 15.9 15.1 17.9 ...
## $ lwr : num 12.4 15.2 13.9 13.1 15.9 ...
## $ upr : num 16.4 19.2 17.8 17 19.8 ...## Melhorar a legenda.
tx <- paste("Umidade de ", levels(da$A), "%", sep="")
lgd <- simpleKey(columns=1, text=tx,
points=FALSE, lines=TRUE,
corner=c(0.03,0.97))
p1 <- xyplot(reng~pota, groups=A, data=da)
p2 <- xyplot(fit~pota, groups=A, data=pred,
type="l", key=lgd,
xlab=expression("Potássio no solo"~(mg~dm^{-3})),
ylab=expression("Rendimento de grãos"~(g~vaso^{-1})),
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.1, fill=1,
prepanel=prepanel.cbH, panel.groups=panel.cbH,
panel=panel.superpose)
p1+as.layer(p2)
## Como calcular os pontos de máximo? São dados por x_max = -0.5*b1/b2.
##-----------------------------------------------------------------------------
## Testar se o rendimento sem potássio (0) é o mesmo para as águas.
X <- LSmatrix(m2, effect="A",
at=list("pota"=0, "I(pota^2)"=0))
rownames(X) <- levels(da$A)
X## (Intercept) blocII blocIII blocIV blocV A50 A62.5 pota I(pota^2) A50:pota A62.5:pota
## 37.5 1 0.2 0.2 0.2 0.2 0 0 0 0 0 0
## 50 1 0.2 0.2 0.2 0.2 1 0 0 0 0 0
## 62.5 1 0.2 0.2 0.2 0.2 0 1 0 0 0 0g <- glht(m2, linfct=X)
confint(g, calpha=univariate_calpha())##
## Simultaneous Confidence Intervals
##
## Fit: lm(formula = reng ~ bloc + A + pota + I(pota^2) + A:pota, data = da)
##
## Quantile = 1.9983
## 95% confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## 37.5 == 0 13.8704 12.2268 15.5141
## 50 == 0 16.6523 15.0086 18.2959
## 62.5 == 0 15.3039 13.6621 16.9458## Contrastes par a par.
Xc <- apc(X)
summary(glht(m2, linfct=Xc))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = reng ~ bloc + A + pota + I(pota^2) + A:pota, data = da)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 37.5:0:0-50:0:0 == 0 -2.782 1.060 -2.625 0.0289 *
## 37.5:0:0-62.5:0:0 == 0 -1.433 1.060 -1.352 0.3720
## 50:0:0-62.5:0:0 == 0 1.348 1.060 1.272 0.4161
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)summary(glht(m2, linfct=Xc), test=Ftest())##
## General Linear Hypotheses
##
## Linear Hypotheses:
## Estimate
## 37.5:0:0-50:0:0 == 0 -2.782
## 37.5:0:0-62.5:0:0 == 0 -1.433
## 50:0:0-62.5:0:0 == 0 1.348
##
## Global Test:
## F DF1 DF2 Pr(>F)
## 1 3.446 2 63 0.03799## O teste simultâneo, intercepto comum.
formula(m1)## reng ~ bloc + A + pota + I(pota^2) + A:pota + A:I(pota^2)m3 <- update(m2, formula=.~.-A)
coef(m3)## (Intercept) blocII blocIII blocIV blocV pota
## 15.8313206978 1.0660000000 -0.8606666667 -1.4433333333 -1.5406099373 0.2022515736
## I(pota^2) pota:A50 pota:A62.5
## -0.0008561129 0.0291812865 0.0742331466anova(m3, m2)## Analysis of Variance Table
##
## Model 1: reng ~ bloc + pota + I(pota^2) + A:pota
## Model 2: reng ~ bloc + A + pota + I(pota^2) + A:pota
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 65 399.38
## 2 63 359.99 2 39.384 3.4462 0.03799 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2)## Analysis of Variance Table
##
## Response: reng
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 4 100.39 25.10 4.392 0.003397 **
## A 2 446.53 223.26 39.072 9.223e-12 ***
## pota 1 1999.89 1999.89 349.988 < 2.2e-16 ***
## I(pota^2) 1 555.07 555.07 97.138 2.200e-14 ***
## A:pota 2 245.17 122.58 21.453 7.841e-08 ***
## Residuals 63 359.99 5.71
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Coeficientes das equações. Mudar a formula do modelo para ter uma
## interpretação de parâmetros por estrato.
## Usar contraste tipo soma para blocos.
formula(m2)## reng ~ bloc + A + pota + I(pota^2) + A:potam3 <- lm(reng~0+A/pota+I(pota^2)+bloc, data=da,
contrasts=list(bloc=contr.sum))
summary(m3)##
## Call:
## lm(formula = reng ~ 0 + A/pota + I(pota^2) + bloc, data = da,
## contrasts = list(bloc = contr.sum))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.7702 -1.4154 -0.0617 1.5246 5.1096
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## A37.5 1.387e+01 8.225e-01 16.864 < 2e-16 ***
## A50 1.665e+01 8.225e-01 20.246 < 2e-16 ***
## A62.5 1.530e+01 8.216e-01 18.627 < 2e-16 ***
## I(pota^2) -8.561e-04 8.602e-05 -9.952 1.51e-14 ***
## bloc1 5.557e-01 5.531e-01 1.005 0.31890
## bloc2 1.622e+00 5.531e-01 2.932 0.00469 **
## bloc3 -3.050e-01 5.531e-01 -0.551 0.58331
## bloc4 -8.876e-01 5.531e-01 -1.605 0.11353
## A37.5:pota 2.129e-01 1.732e-02 12.293 < 2e-16 ***
## A50:pota 2.210e-01 1.732e-02 12.757 < 2e-16 ***
## A62.5:pota 2.763e-01 1.748e-02 15.810 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.39 on 63 degrees of freedom
## Multiple R-squared: 0.9927, Adjusted R-squared: 0.9914
## F-statistic: 778.4 on 11 and 63 DF, p-value: < 2.2e-16## Prova de que são o mesmo modelo.
deviance(m3); deviance(m2)## [1] 359.9933## [1] 359.9933## O R² é o quadrado da correlação entre observado e predito.
cor(da$reng, fitted(m2))^2## [1] 0.9028891## R² separado por água.
ddply(data.frame(o=da$reng, f=fitted(m2)), .(da$A),
summarise, R2=cor(o, f)^2)## da$A R2
## 1 37.5 0.7851939
## 2 50 0.8522193
## 3 62.5 0.9457677##-----------------------------------------------------------------------------
## Os pontos de máximo.
## Buscas essas ocorrências de texto.
oc <- paste0("^A", levels(da$A), ":")
## Estimativas e suas posições de ordem de efeito.
c3 <- coef(m3)
a <- m3$assign
cbind(c3, a) ## c3 a
## A37.5 13.8704427563 1
## A50 16.6522875839 1
## A62.5 15.3039401494 1
## I(pota^2) -0.0008561256 2
## bloc1 0.5556840320 3
## bloc2 1.6216840320 3
## bloc3 -0.3049826347 3
## bloc4 -0.8876493013 3
## A37.5:pota 0.2129359777 4
## A50:pota 0.2209687363 4
## A62.5:pota 0.2762725227 4## Estimativas do ponto de máximo.
ptmax <- -0.5*c3[a==4]/c3[a==2]
xyplot(fit~pota, groups=A, data=pred,
type="l", key=lgd,
xlab=expression("Potássio no solo"~(mg~dm^{-3})),
ylab=expression("Rendimento de grãos"~(g~vaso^{-1})),
ly=pred$lwr, uy=pred$upr, cty="bands", alpha=0.1, fill=1,
prepanel=prepanel.cbH,
panel.groups=panel.cbH,
panel=panel.superpose)+
layer(
panel.abline(v=ptmax, lty=2))
##-----------------------------------------------------------------------------
## E se ajustar um modelo linear-platô? E se for um quadrático-platô?Exemplo 3 (praticar)
##-----------------------------------------------------------------------------
## Dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/algodão.txt",
header=TRUE, sep="\t", encoding="latin1")
da$desf <- da$desf/100
levels(da$estag) <- c("Vegetativo", "Botão floral",
"Florescimento", "Maçã", "Capulho")
str(da)## 'data.frame': 125 obs. of 8 variables:
## $ estag : Factor w/ 5 levels "Vegetativo","Botão floral",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ desf : num 0 0 0 0 0 0.25 0.25 0.25 0.25 0.25 ...
## $ rep : int 1 2 3 4 5 1 2 3 4 5 ...
## $ pcapu : num 33.2 28.7 31.5 28.9 36.4 ...
## $ nnos : int 32 29 29 27 28 30 33 23 29 26 ...
## $ alt : num 151 146 140 155 150 ...
## $ ncapu : int 10 9 8 8 10 11 9 10 10 10 ...
## $ nestru: int 10 9 8 9 10 11 9 11 10 10 ...xyplot(pcapu~desf|estag, data=da, type=c("p","smooth"),
xlab="Nível de desfolha artificial",
ylab="Peso de capulhos produzidos (g)")
Análise de covariância
##-----------------------------------------------------------------------------
## Lendo arquivos de dados.
## Url de um arquivo com dados.
url <- "http://www.leg.ufpr.br/~walmes/data/ancova.txt"
## Importa os dados a partir do arquivo na web.
da <- read.table(file=url, header=TRUE, sep="\t")
## Trunca os nomes para 4 digitos.
names(da) <- substr(names(da), 1, 4)
## Mostra a estrutura do objeto.
str(da)## 'data.frame': 72 obs. of 8 variables:
## $ ener: Factor w/ 3 levels "alto","baixo",..: 2 2 2 2 2 2 2 2 3 3 ...
## $ sexo: Factor w/ 3 levels "F","MC","MI": 2 2 2 2 2 2 2 2 2 2 ...
## $ rep : int 1 2 3 4 5 6 7 8 1 2 ...
## $ pi : num 93.3 94 91.6 89.3 91.5 93 88.8 93.7 93.9 88.5 ...
## $ id : int 134 137 133 133 144 147 138 142 134 137 ...
## $ pviv: num 127 125 126 119 122 ...
## $ rend: num 82.3 83.1 80.7 82.7 82.9 ...
## $ peso: num 127 125 126 119 122 ...## Dados de experimento com nutrição de suínos. Animais foram pesados
## antes do experimento e tinham idade conhecida. Essas variáveis
## contínuas foram usadas para explicar/corrigir parte da variação
## presente e melhor comparar os níveis de energia na ração fornecidos.\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Normal}(\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \mu+\alpha_i+\tau_j+\gamma_{ij} \end{aligned} \]
em que \(\alpha\) é o efeito dos níveis de sexo \(i\), \(\tau_j\) o efeito dos níveis de energia (Kcal) na dieta \(j\), \(gamma_{ij}\) representa a interação . \(\mu\) é a média de \(Y\) na ausência do efeito dos termos mencioandos. \(\mu_{ij}\) é a média para as observações na combinação \(ij\) e \(\sigma^2\) é a variância das observações ao redor desse valor médio. No modelo de análise de covariância, os valores de peso inicial e idade serão usados para explicar parte da variação ao final do experimento. Sendo assim, o modelo considerando o efeito dessas variáveis é representado por
\[ \begin{aligned} Y|\text{fontes de variação} &\,\sim \text{Normal}(\mu+\beta_1 x+\beta_2 z+\mu_{ij}, \sigma^2) \newline \mu_{ij} &= \alpha_i+\tau_j+\gamma_{ij} \end{aligned} \]
em que \(\beta_1\) representa o efeito da covariável peso (\(x\)) expresso por um coeficiente angular e \(\beta_2\) o mesmo para o efeito da covariável idade inicial (\(z\)).
##-----------------------------------------------------------------------------
## Especificação dos modelos.
## Só com as fontes de variação de interesse.
m0 <- lm(peso~sexo*ener, data=da)
anova(m0)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## sexo 2 199.81 99.903 3.6153 0.03263 *
## ener 2 55.84 27.921 1.0104 0.36990
## sexo:ener 4 35.38 8.846 0.3201 0.86348
## Residuals 63 1740.92 27.634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Análise de variância (em experimentos não ortogonais a ordem dos
## termos é importante!).
m1 <- lm(peso~pi+id+sexo*ener, data=da)
anova(m1)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## pi 1 595.58 595.58 32.1128 4.209e-07 ***
## id 1 47.23 47.23 2.5464 0.11571
## sexo 2 156.97 78.48 4.2318 0.01901 *
## ener 2 66.46 33.23 1.7918 0.17531
## sexo:ener 4 34.38 8.59 0.4634 0.76228
## Residuals 61 1131.33 18.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Troca ordem de entrada apenas por curiosidade.
m1 <- lm(peso~id+pi+sexo*ener, data=da)
anova(m1)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## id 1 5.17 5.17 0.2789 0.59937
## pi 1 637.63 637.63 34.3803 1.98e-07 ***
## sexo 2 156.97 78.48 4.2318 0.01901 *
## ener 2 66.46 33.23 1.7918 0.17531
## sexo:ener 4 34.38 8.59 0.4634 0.76228
## Residuals 61 1131.33 18.55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Testa o poder de explicação dos "blocos" contínuos, termos extras.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: peso ~ sexo * ener
## Model 2: peso ~ id + pi + sexo * ener
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 63 1740.9
## 2 61 1131.3 2 609.59 16.434 1.953e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m1); layout(1)
## Medidas de influência para as observações.
im <- influence.measures(m1)
summary(im)## Potentially influential observations of
## lm(formula = peso ~ id + pi + sexo * ener, data = da) :
##
## dfb.1_ dfb.id dfb.pi dfb.sxMC dfb.sxMI dfb.enrb dfb.enrm dfb.sxMC:nrb dfb.sxMI:nrb
## 25 -0.60 -0.42 0.95 -0.04 0.00 -0.02 0.02 0.05 0.68
## 37 -0.20 0.96 -0.43 0.07 0.00 0.01 -0.04 -0.04 0.09
## 46 0.40 -0.57 -0.05 -0.04 -0.82 0.00 0.02 0.01 0.53
## 48 0.28 0.30 -0.52 0.03 0.81 0.01 -0.01 -0.03 -0.54
## dfb.sxMC:nrm dfb.sxMI:nrm dffit cov.r cook.d hat
## 25 -0.01 -0.06 1.76_* 0.15_* 0.23 0.18
## 37 0.12 -0.43 -1.48_* 0.44_* 0.18 0.23
## 46 -0.09 0.51 -1.30_* 0.30_* 0.14 0.16
## 48 0.02 -0.53 1.27 0.31_* 0.13 0.15r <- which(im$is.inf[,"dffit"])
m2 <- lm(peso~id+pi+sexo*ener, data=da[-r,])
anova(m2)## Analysis of Variance Table
##
## Response: peso
## Df Sum Sq Mean Sq F value Pr(>F)
## id 1 0.40 0.40 0.0341 0.854171
## pi 1 459.25 459.25 39.3607 4.8e-08 ***
## sexo 2 177.89 88.94 7.6230 0.001150 **
## ener 2 144.97 72.49 6.2125 0.003592 **
## sexo:ener 4 66.83 16.71 1.4319 0.234964
## Residuals 58 676.73 11.67
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1par(mfrow=c(2,2)); plot(m2); layout(1)
##-----------------------------------------------------------------------------
## Comparações múltiplas para o desdobramento.
LSmeans(m2, effect=c("sexo","ener"))## estimate se df t.stat p.value sexo ener id pi
## 1 126.9699 1.211308 58 104.82046 8.066598e-68 F alto 138.0145 92.0087
## 2 124.9146 1.209095 58 103.31250 1.861161e-67 MC alto 138.0145 92.0087
## 3 129.7600 1.291111 58 100.50258 9.132251e-67 MI alto 138.0145 92.0087
## 4 122.7648 1.211616 58 101.32318 5.713283e-67 F baixo 138.0145 92.0087
## 5 124.0100 1.208737 58 102.59473 2.782637e-67 MC baixo 138.0145 92.0087
## 6 124.5507 1.293884 58 96.26108 1.097349e-65 MI baixo 138.0145 92.0087
## 7 125.6840 1.219403 58 103.07010 2.131268e-67 F medio 138.0145 92.0087
## 8 123.9038 1.251383 58 99.01350 2.159913e-66 MC medio 138.0145 92.0087
## 9 130.0355 1.291530 58 100.68331 8.233295e-67 MI medio 138.0145 92.0087LSmeans(m2, effect=c("sexo"))## estimate se df t.stat p.value sexo id pi
## 1 125.1396 0.7071868 58 176.9540 5.740376e-81 F 138.0145 92.0087
## 2 124.2761 0.7048357 58 176.3193 7.067888e-81 MC 138.0145 92.0087
## 3 128.1154 0.7456101 58 171.8262 3.149698e-80 MI 138.0145 92.0087LSmeans(m2, effect=c("ener"))## estimate se df t.stat p.value ener id pi
## 1 127.2148 0.7141424 58 178.1365 3.903753e-81 alto 138.0145 92.0087
## 2 123.7752 0.7138463 58 173.3919 1.863060e-80 baixo 138.0145 92.0087
## 3 126.5411 0.7149777 58 176.9861 5.680494e-81 medio 138.0145 92.0087Xc <- LSmatrix(m2, effect=c("sexo"))
rownames(X) <- levels(da$sexo)
ps <- apmc(X, model=m2, focus="sexo", test="fdr")
ps## sexo estimate lwr upr cld
## 1 F 40.85947 4.41751429 77.30142 a
## 2 MC 36.65442 0.09577609 73.21306 a
## 3 MI 39.57363 2.99722792 76.15004 aX <- LSmatrix(m2, effect=c("ener"))
rownames(X) <- levels(da$ener)
pe <- apmc(X, model=m2, focus="ener", test="fdr")
pe$ener <- factor(pe$ener, levels=c("alto","medio","baixo"),
labels=c("Alto","Médio","Baixo"))
pe <- arrange(pe, ener)
pe## ener estimate lwr upr cld
## 1 Alto 127.2148 125.7853 128.6443 a
## 2 Médio 126.5411 125.1099 127.9723 a
## 3 Baixo 123.7752 122.3462 125.2041 b## names(ps)[1] <- names(pe)[1] <- "nivel"
## p0 <- cbind(rbind(ps, pe), fator=rep(c("sexo","ener"), each=3))
## dput(levels(p0$nivel))
## str(p0)##-----------------------------------------------------------------------------
## Representação gráfica dos resultados.
xlab <- expression("Peso aos 28 dias"~(kg))
sub <- list(
"Médias seguidas de mesma letra indicam diferença nula à 5%.",
font=1, cex=0.8)
p1 <-
segplot(sexo~lwr+upr, centers=estimate,
data=ps, draw=FALSE,
ylab="Nível de sexo", xlab=xlab, sub=sub,
letras=sprintf("%0.2f %s", ps$estimate, ps$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(y=as.numeric(z)[subscripts],
x=centers[subscripts],
label=letras[subscripts], pos=3)
})
p2 <-
segplot(ener~lwr+upr, centers=estimate,
data=pe, draw=FALSE,
ylab="Nível de energia da dieta", xlab=xlab, sub=sub,
letras=sprintf("%0.2f %s", pe$estimate, pe$cld),
panel=function(x, y, z, subscripts, centers, letras, ...){
panel.segplot(x=x, y=y, z=z, centers=centers,
subscripts=subscripts, ...)
panel.text(y=as.numeric(z)[subscripts],
x=centers[subscripts],
label=letras[subscripts], pos=3)
})
## Gráfico final.
plot(p1, split=c(1,1,2,1), more=TRUE)
plot(p2, split=c(2,1,2,1), more=FALSE)
##-----------------------------------------------------------------------------
## Como fica o resultado sem usar as covariáveis?
mean(da$pi); mean(da$id)## [1] 92.11667## [1] 137.8889X <- LSmatrix(m0, effect=c("sexo"), at=list(pi=92, id=138))
rownames(X) <- levels(da$sexo)
ps0 <- apmc(X, model=m0, focus="sexo", test="fdr")
ps0## sexo estimate lwr upr cld
## 1 F 125.1167 122.9724 127.2610 ab
## 2 MC 124.3250 122.1807 126.4693 b
## 3 MI 128.1875 126.0432 130.3318 aX <- LSmatrix(m1, effect=c("sexo"), at=list(pi=92, id=138))
rownames(X) <- levels(da$sexo)
ps1 <- apmc(X, model=m1, focus="sexo", test="fdr")
ps1## sexo estimate lwr upr cld
## 1 F 125.1306 123.3506 126.9106 ab
## 2 MC 124.2653 122.4914 126.0392 b
## 3 MI 127.7405 125.9733 129.5076 aX <- LSmatrix(m2, effect=c("sexo"), at=list(pi=92, id=138))
rownames(X) <- levels(da$sexo)
ps2 <- apmc(X, model=m2, focus="sexo", test="fdr")
ps2## sexo estimate lwr upr cld
## 1 F 125.1316 123.7156 126.5476 b
## 2 MC 124.2682 122.8576 125.6787 b
## 3 MI 128.1074 126.6150 129.5999 aps0$model <- "0"; ps1$model <- "1"; ps2$model <- "2"
ps <- rbind(ps0, ps1, ps2)
ps$model <- factor(ps$model)
str(ps)## 'data.frame': 9 obs. of 6 variables:
## $ sexo : Factor w/ 3 levels "F","MC","MI": 1 2 3 1 2 3 1 2 3
## $ estimate: num 125 124 128 125 124 ...
## $ lwr : num 123 122 126 123 122 ...
## $ upr : num 127 126 130 127 126 ...
## $ cld : chr "ab" "b" "a" "ab" ...
## $ model : Factor w/ 3 levels "0","1","2": 1 1 1 2 2 2 3 3 3l <- c("Fatorial", "Ancova", "Limpo")
key <- list(title="Modelo", cex.title=1.1,
corner=c(0.05,0.95), type="o", divide=1,
text=list(l), lines=list(pch=1:length(l)))
segplot(sexo~lwr+upr, centers=estimate,
data=ps, groups=model, draw=FALSE,
ylab="Nível de sexo", xlab=xlab, sub=sub, key=key,
panel=panel.segplotBy, f=0.075, pch=as.integer(ps$model))
Experimento fatorial triplo com regressão
http://www.leg.ufpr.br/~walmes/analises/MESerafim/maracuja/
##-----------------------------------------------------------------------------
## Ler tabela de dados.
da <- read.table("http://www.leg.ufpr.br/~walmes/data/maracuja_planta.txt",
header=TRUE, sep="\t")
da$bloc <- factor(da$bloc)
str(da)## 'data.frame': 168 obs. of 9 variables:
## $ agreg: Factor w/ 2 levels "[0,2]","[4,10]": 2 2 2 2 2 2 2 2 2 2 ...
## $ fam : Factor w/ 2 levels "F29","F48": 1 1 1 1 1 1 1 1 1 1 ...
## $ cinz : num 0 0 0 0 0 0 1.5 1.5 1.5 1.5 ...
## $ bloc : Factor w/ 3 levels "1","2","3": 1 1 2 2 3 3 1 1 2 2 ...
## $ rept : int 1 2 1 2 1 2 1 2 1 2 ...
## $ mfpa : num 43.5 43.7 31.8 23.4 36.4 ...
## $ mspa : num 11.1 11.13 8.24 5.39 9.04 ...
## $ ds : num 1.24 1.25 1.23 1.23 1.27 1.18 1.27 1.28 1.01 1.19 ...
## $ cav : num 6470 7260 9446 9046 8468 ...u <- unique(da$cinz)
cbind(u, u/1.5, sqrt(u/1.5), log(u/1.5, base=2))## u
## [1,] 0.0 0 0.000000 -Inf
## [2,] 1.5 1 1.000000 0
## [3,] 3.0 2 1.414214 1
## [4,] 6.0 4 2.000000 2
## [5,] 12.0 8 2.828427 3
## [6,] 24.0 16 4.000000 4
## [7,] 48.0 32 5.656854 5## Doses de cinza em uma escala com distância mais regular entre níveis.
da$cin <- sqrt(da$cinz/1.5)
##-----------------------------------------------------------------------------
## Ver.
p1 <- xyplot(mfpa~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
p2 <- xyplot(mspa~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
p3 <- xyplot(ds~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
p4 <- xyplot(cav~cin|fam, groups=agreg, data=da, type=c("p","smooth"),
auto.key=list(columns=2))
grid.arrange(p1, p2, nrow=2)
grid.arrange(p3, p4, nrow=2)
##-----------------------------------------------------------------------------
## Especificação e ajuste dos modelos em batelada.
resp <- c("mfpa", "mspa", "ds", "cav")
form0 <- lapply(paste(resp, "~bloc+fam*agreg*(cinz+I(cinz^2))"),
as.formula)
names(form0) <- resp
## Ajustes.
m0 <- lapply(form0, lm, data=da)
##-----------------------------------------------------------------------------
## Quadros de anova.
lapply(m0, anova)## $mfpa
## Analysis of Variance Table
##
## Response: mfpa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 2103.5 1051.8 6.5549 0.001853 **
## fam 1 88.7 88.7 0.5525 0.458424
## agreg 1 4129.9 4129.9 25.7388 1.111e-06 ***
## cinz 1 1848.3 1848.3 11.5191 0.000876 ***
## I(cinz^2) 1 426.3 426.3 2.6569 0.105148
## fam:agreg 1 257.8 257.8 1.6065 0.206896
## fam:cinz 1 431.2 431.2 2.6872 0.103198
## fam:I(cinz^2) 1 559.1 559.1 3.4847 0.063841 .
## agreg:cinz 1 839.9 839.9 5.2348 0.023500 *
## agreg:I(cinz^2) 1 26.5 26.5 0.1654 0.684839
## fam:agreg:cinz 1 252.3 252.3 1.5726 0.211733
## fam:agreg:I(cinz^2) 1 122.8 122.8 0.7652 0.383074
## Residuals 154 24710.0 160.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $mspa
## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 86.24 43.121 5.2662 0.006134 **
## fam 1 0.32 0.315 0.0385 0.844648
## agreg 1 80.26 80.261 9.8020 0.002086 **
## cinz 1 61.76 61.758 7.5423 0.006745 **
## I(cinz^2) 1 4.97 4.970 0.6069 0.437143
## fam:agreg 1 1.24 1.245 0.1520 0.697172
## fam:cinz 1 20.51 20.512 2.5051 0.115529
## fam:I(cinz^2) 1 21.27 21.272 2.5979 0.109053
## agreg:cinz 1 63.03 63.027 7.6973 0.006215 **
## agreg:I(cinz^2) 1 13.61 13.609 1.6620 0.199267
## fam:agreg:cinz 1 11.24 11.240 1.3728 0.243147
## fam:agreg:I(cinz^2) 1 7.39 7.387 0.9021 0.343697
## Residuals 154 1260.98 8.188
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $ds
## Analysis of Variance Table
##
## Response: ds
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 0.04405 0.02203 5.0632 0.0074192 **
## fam 1 0.05306 0.05306 12.1962 0.0006252 ***
## agreg 1 0.02394 0.02394 5.5037 0.0202505 *
## cinz 1 0.81520 0.81520 187.3804 < 2.2e-16 ***
## I(cinz^2) 1 0.08952 0.08952 20.5770 1.146e-05 ***
## fam:agreg 1 0.01259 0.01259 2.8940 0.0909284 .
## fam:cinz 1 0.03576 0.03576 8.2192 0.0047257 **
## fam:I(cinz^2) 1 0.03236 0.03236 7.4382 0.0071268 **
## agreg:cinz 1 0.02580 0.02580 5.9314 0.0160155 *
## agreg:I(cinz^2) 1 0.01699 0.01699 3.9060 0.0499013 *
## fam:agreg:cinz 1 0.00661 0.00661 1.5188 0.2196819
## fam:agreg:I(cinz^2) 1 0.00008 0.00008 0.0187 0.8914281
## Residuals 154 0.66998 0.00435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $cav
## Analysis of Variance Table
##
## Response: cav
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 4704295 2352148 1.5448 0.2166456
## fam 1 12968149 12968149 8.5169 0.0040462 **
## agreg 1 270606640 270606640 177.7222 < 2.2e-16 ***
## cinz 1 550533 550533 0.3616 0.5485225
## I(cinz^2) 1 2360352 2360352 1.5502 0.2150012
## fam:agreg 1 4240215 4240215 2.7848 0.0971955 .
## fam:cinz 1 581539 581539 0.3819 0.5374858
## fam:I(cinz^2) 1 14011881 14011881 9.2024 0.0028376 **
## agreg:cinz 1 2821636 2821636 1.8531 0.1754091
## agreg:I(cinz^2) 1 13641 13641 0.0090 0.9247148
## fam:agreg:cinz 1 19103685 19103685 12.5464 0.0005256 ***
## fam:agreg:I(cinz^2) 1 291877 291877 0.1917 0.6621261
## Residuals 154 234486288 1522638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##-----------------------------------------------------------------------------
## Quadro de estimativas dos efeitos.
lapply(m0, summary)## $mfpa
##
## Call:
## FUN(formula = X[[1L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.074 -8.108 1.068 7.509 32.862
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.947583 3.468766 11.516 < 2e-16 ***
## bloc2 -2.623214 2.393848 -1.096 0.274871
## bloc3 -8.465893 2.393848 -3.537 0.000536 ***
## famF48 8.487977 4.499371 1.886 0.061113 .
## agreg[4,10] -0.716554 4.499371 -0.159 0.873676
## cinz 1.212067 0.464685 2.608 0.009993 **
## I(cinz^2) -0.018792 0.009479 -1.982 0.049212 *
## famF48:agreg[4,10] -12.400527 6.363071 -1.949 0.053133 .
## famF48:cinz -1.182186 0.657164 -1.799 0.073989 .
## famF48:I(cinz^2) 0.025987 0.013406 1.939 0.054389 .
## agreg[4,10]:cinz -0.642157 0.657164 -0.977 0.330020
## agreg[4,10]:I(cinz^2) 0.004437 0.013406 0.331 0.741092
## famF48:agreg[4,10]:cinz 1.090510 0.929370 1.173 0.242452
## famF48:agreg[4,10]:I(cinz^2) -0.016584 0.018958 -0.875 0.383074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.67 on 154 degrees of freedom
## Multiple R-squared: 0.3097, Adjusted R-squared: 0.2514
## F-statistic: 5.315 on 13 and 154 DF, p-value: 7.802e-08
##
##
## $mspa
##
## Call:
## FUN(formula = X[[2L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3162 -1.8967 0.2143 2.0662 6.9304
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.3129949 0.7835985 11.885 < 2e-16 ***
## bloc2 -0.2450000 0.5407732 -0.453 0.65115
## bloc3 -1.6275000 0.5407732 -3.010 0.00306 **
## famF48 1.0284923 1.0164134 1.012 0.31318
## agreg[4,10] 0.1132211 1.0164134 0.111 0.91145
## cinz 0.1743295 0.1049729 1.661 0.09881 .
## I(cinz^2) -0.0021965 0.0021414 -1.026 0.30662
## famF48:agreg[4,10] -2.0292102 1.4374256 -1.412 0.16006
## famF48:cinz -0.2481261 0.1484541 -1.671 0.09667 .
## famF48:I(cinz^2) 0.0054853 0.0030283 1.811 0.07204 .
## agreg[4,10]:cinz -0.0744156 0.1484541 -0.501 0.61690
## agreg[4,10]:I(cinz^2) -0.0007267 0.0030283 -0.240 0.81067
## famF48:agreg[4,10]:cinz 0.2570117 0.2099458 1.224 0.22275
## famF48:agreg[4,10]:I(cinz^2) -0.0040678 0.0042827 -0.950 0.34370
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.862 on 154 degrees of freedom
## Multiple R-squared: 0.2277, Adjusted R-squared: 0.1625
## F-statistic: 3.493 on 13 and 154 DF, p-value: 8.922e-05
##
##
## $ds
##
## Call:
## FUN(formula = X[[3L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.271875 -0.030538 0.003811 0.033212 0.179064
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.264e+00 1.806e-02 69.958 < 2e-16 ***
## bloc2 -2.099e-02 1.246e-02 -1.684 0.09417 .
## bloc3 1.865e-02 1.246e-02 1.496 0.13666
## famF48 4.518e-02 2.343e-02 1.929 0.05562 .
## agreg[4,10] -3.503e-02 2.343e-02 -1.495 0.13695
## cinz -8.080e-03 2.420e-03 -3.339 0.00105 **
## I(cinz^2) 9.004e-05 4.936e-05 1.824 0.07005 .
## famF48:agreg[4,10] -5.847e-02 3.313e-02 -1.765 0.07960 .
## famF48:cinz -9.292e-03 3.422e-03 -2.715 0.00738 **
## famF48:I(cinz^2) 1.414e-04 6.980e-05 2.025 0.04458 *
## agreg[4,10]:cinz 5.060e-03 3.422e-03 1.479 0.14130
## agreg[4,10]:I(cinz^2) -9.080e-05 6.980e-05 -1.301 0.19527
## famF48:agreg[4,10]:cinz 2.205e-03 4.839e-03 0.456 0.64933
## famF48:agreg[4,10]:I(cinz^2) -1.350e-05 9.872e-05 -0.137 0.89143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06596 on 154 degrees of freedom
## Multiple R-squared: 0.6331, Adjusted R-squared: 0.6021
## F-statistic: 20.44 on 13 and 154 DF, p-value: < 2.2e-16
##
##
## $cav
##
## Call:
## FUN(formula = X[[4L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2658.0 -767.1 -123.3 800.8 3245.6
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6018.4385 337.9074 17.811 < 2e-16 ***
## bloc2 -10.1964 233.1951 -0.044 0.96518
## bloc3 349.7679 233.1951 1.500 0.13569
## famF48 -1070.2391 438.3031 -2.442 0.01575 *
## agreg[4,10] 1370.2799 438.3031 3.126 0.00212 **
## cinz -130.1812 45.2670 -2.876 0.00460 **
## I(cinz^2) 2.1339 0.9234 2.311 0.02217 *
## famF48:agreg[4,10] 1934.0225 619.8542 3.120 0.00216 **
## famF48:cinz 201.0797 64.0171 3.141 0.00202 **
## famF48:I(cinz^2) -3.2055 1.3059 -2.455 0.01522 *
## agreg[4,10]:cinz 73.2835 64.0171 1.145 0.25409
## agreg[4,10]:I(cinz^2) -0.3169 1.3059 -0.243 0.80859
## famF48:agreg[4,10]:cinz -122.4676 90.5339 -1.353 0.17813
## famF48:agreg[4,10]:I(cinz^2) 0.8086 1.8468 0.438 0.66213
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1234 on 154 degrees of freedom
## Multiple R-squared: 0.5863, Adjusted R-squared: 0.5513
## F-statistic: 16.79 on 13 and 154 DF, p-value: < 2.2e-16##-----------------------------------------------------------------------------
## Avaliação dos pressupostos.
par(mfrow=c(2,2)); plot(m0[["mfpa"]]); layout(1) ## ok!
par(mfrow=c(2,2)); plot(m0[["mspa"]]); layout(1) ## ok!
par(mfrow=c(2,2)); plot(m0[["ds"]]); layout(1) ## bom.
par(mfrow=c(2,2)); plot(m0[["cav"]]); layout(1) ## ok!
## par(mfrow=c(2,2))
## lapply(m0, plot, which=1)
## layout(1)
## par(mfrow=c(2,2))
## lapply(m0, plot, which=2)
## layout(1)
## par(mfrow=c(2,2))
## lapply(m0, plot, which=3)
## layout(1)
##-----------------------------------------------------------------------------
## Modelos após abandono de termos não significativos.
form1 <- list(mfpa=mfpa~bloc+agreg*cinz,
mspa=mspa~bloc+agreg*cinz,
ds=ds~bloc+(fam+agreg+cinz)^3+(fam+agreg)*I(cinz^2),
cav=cav~bloc+fam*agreg*(cinz+I(cinz^2)))
m1 <- lapply(form1, lm, data=da, contrast=list(bloc=contr.sum))
lapply(m1, anova)## $mfpa
## Analysis of Variance Table
##
## Response: mfpa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 2103.5 1051.8 6.3401 0.002234 **
## agreg 1 4129.9 4129.9 24.8950 1.548e-06 ***
## cinz 1 1848.3 1848.3 11.1414 0.001047 **
## agreg:cinz 1 839.9 839.9 5.0632 0.025786 *
## Residuals 162 26874.6 165.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $mspa
## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 86.24 43.121 5.2072 0.006431 **
## agreg 1 80.26 80.261 9.6921 0.002188 **
## cinz 1 61.76 61.758 7.4577 0.007016 **
## agreg:cinz 1 63.03 63.027 7.6110 0.006469 **
## Residuals 162 1341.53 8.281
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $ds
## Analysis of Variance Table
##
## Response: ds
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 0.04405 0.02203 5.0955 0.0071909 **
## fam 1 0.05306 0.05306 12.2739 0.0006005 ***
## agreg 1 0.02394 0.02394 5.5388 0.0198546 *
## cinz 1 0.81520 0.81520 188.5742 < 2.2e-16 ***
## I(cinz^2) 1 0.08952 0.08952 20.7081 1.074e-05 ***
## fam:agreg 1 0.01259 0.01259 2.9124 0.0899018 .
## fam:cinz 1 0.03576 0.03576 8.2716 0.0045943 **
## agreg:cinz 1 0.02580 0.02580 5.9691 0.0156812 *
## fam:I(cinz^2) 1 0.03236 0.03236 7.4856 0.0069454 **
## agreg:I(cinz^2) 1 0.01699 0.01699 3.9308 0.0491754 *
## fam:agreg:cinz 1 0.00661 0.00661 1.5285 0.2182139
## Residuals 155 0.67006 0.00432
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $cav
## Analysis of Variance Table
##
## Response: cav
## Df Sum Sq Mean Sq F value Pr(>F)
## bloc 2 4704295 2352148 1.5448 0.2166456
## fam 1 12968149 12968149 8.5169 0.0040462 **
## agreg 1 270606640 270606640 177.7222 < 2.2e-16 ***
## cinz 1 550533 550533 0.3616 0.5485225
## I(cinz^2) 1 2360352 2360352 1.5502 0.2150012
## fam:agreg 1 4240215 4240215 2.7848 0.0971955 .
## fam:cinz 1 581539 581539 0.3819 0.5374858
## fam:I(cinz^2) 1 14011881 14011881 9.2024 0.0028376 **
## agreg:cinz 1 2821636 2821636 1.8531 0.1754091
## agreg:I(cinz^2) 1 13641 13641 0.0090 0.9247148
## fam:agreg:cinz 1 19103685 19103685 12.5464 0.0005256 ***
## fam:agreg:I(cinz^2) 1 291877 291877 0.1917 0.6621261
## Residuals 154 234486288 1522638
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Anova entre modelos sequenciais.
sapply(names(m1), simplify=FALSE,
function(i){
anova(m0[[i]],m1[[i]])
})## $mfpa
## Analysis of Variance Table
##
## Model 1: mfpa ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: mfpa ~ bloc + agreg * cinz
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 24710
## 2 162 26875 -8 -2164.7 1.6864 0.1058
##
## $mspa
## Analysis of Variance Table
##
## Model 1: mspa ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: mspa ~ bloc + agreg * cinz
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 1261.0
## 2 162 1341.5 -8 -80.55 1.2297 0.2852
##
## $ds
## Analysis of Variance Table
##
## Model 1: ds ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: ds ~ bloc + (fam + agreg + cinz)^3 + (fam + agreg) * I(cinz^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 0.66998
## 2 155 0.67006 -1 -8.1325e-05 0.0187 0.8914
##
## $cav
## Analysis of Variance Table
##
## Model 1: cav ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Model 2: cav ~ bloc + fam * agreg * (cinz + I(cinz^2))
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 154 234486288
## 2 154 234486288 0 2.9802e-08##-----------------------------------------------------------------------------
##
summary(m0[["mfpa"]])##
## Call:
## FUN(formula = X[[1L]], data = ..1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.074 -8.108 1.068 7.509 32.862
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.947583 3.468766 11.516 < 2e-16 ***
## bloc2 -2.623214 2.393848 -1.096 0.274871
## bloc3 -8.465893 2.393848 -3.537 0.000536 ***
## famF48 8.487977 4.499371 1.886 0.061113 .
## agreg[4,10] -0.716554 4.499371 -0.159 0.873676
## cinz 1.212067 0.464685 2.608 0.009993 **
## I(cinz^2) -0.018792 0.009479 -1.982 0.049212 *
## famF48:agreg[4,10] -12.400527 6.363071 -1.949 0.053133 .
## famF48:cinz -1.182186 0.657164 -1.799 0.073989 .
## famF48:I(cinz^2) 0.025987 0.013406 1.939 0.054389 .
## agreg[4,10]:cinz -0.642157 0.657164 -0.977 0.330020
## agreg[4,10]:I(cinz^2) 0.004437 0.013406 0.331 0.741092
## famF48:agreg[4,10]:cinz 1.090510 0.929370 1.173 0.242452
## famF48:agreg[4,10]:I(cinz^2) -0.016584 0.018958 -0.875 0.383074
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 12.67 on 154 degrees of freedom
## Multiple R-squared: 0.3097, Adjusted R-squared: 0.2514
## F-statistic: 5.315 on 13 and 154 DF, p-value: 7.802e-08##-----------------------------------------------------------------------------
## Valores preditos.
gridlist <- list(bloc="1",
fam=levels(da$fam),
agreg=levels(da$agreg),
cin=seq(0,6,l=100),
cinz=seq(0,50,l=100))
m1 <- lapply(m1,
function(i){
predvars <- attr(terms(i), "term.labels")
pred <- do.call(
expand.grid,
gridlist[predvars[!sapply(gridlist[predvars],
is.null)]])
i$newdata <- pred
i
})
mypredict <- function(m){
cbind(m$newdata,
predict(m, newdata=m$newdata,
interval="confidence")-
coef(m)["bloc1"])
}
all.pred <- lapply(m1, mypredict)
## str(all.pred)
## lapply(all.pred, names)
##-----------------------------------------------------------------------------
## Gráficos.
xlab <- expression(Cinza~(t~ha^{-1}))
ylab <- list(expression(Massa~fresca~de~parte~aérea~(g)),
expression(Massa~seca~de~parte~aérea~(g)),
expression(Densidade~do~solo~(Mg~t^{-1})),
expression(Água~consumida~no~ciclo~(mL)))
names(ylab) <- names(m1)
##-----------------------------------------------------------------------------
## Mfpa.
p1 <-
xyplot(mfpa~cinz, groups=agreg, data=da,
ylab=ylab[["mfpa"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["mfpa"]],
xyplot(fit~cinz, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Mspa.
p2 <-
xyplot(mspa~cinz, groups=agreg, data=da,
ylab=ylab[["mspa"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["mspa"]],
xyplot(fit~cinz, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Ds.
p3 <-
xyplot(ds~cinz|fam, groups=agreg, data=da, layout=c(NA,1),
ylab=ylab[["ds"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["ds"]],
xyplot(fit~cinz|fam, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Cav.
p4 <-
xyplot(cav~cinz|fam, groups=agreg, data=da, layout=c(NA,1),
ylab=ylab[["cav"]], xlab=xlab,
strip=strip.custom(bg="gray90"),
auto.key=list(columns=2, points=FALSE, lines=TRUE,
title="Classe de agregado (mm)", cex.title=1.2))+
as.layer(with(all.pred[["cav"]],
xyplot(fit~cinz|fam, groups=agreg, type="l",
ly=lwr, uy=upr, cty="bands", fill=1, alpha=0.2,
prepanel=prepanel.cbH,
panel=panel.superpose,
panel.groups=panel.cbH)))
##-----------------------------------------------------------------------------
## Arranjo 1.
grid.arrange(p1, p2, ncol=2)
grid.text(label="A", x=0.025, y=0.975)
grid.text(label="B", x=0.5+0.025, y=0.975)
##-----------------------------------------------------------------------------
## Arranjo 2.
grid.arrange(p3, p4, ncol=2)
grid.text(label="A", x=0.025, y=0.975)
grid.text(label="B", x=0.5+0.025, y=0.975)
##-----------------------------------------------------------------------------
## Informações da sessão.
Sys.time()## [1] "2014-12-11 20:59:05 BRST"sessionInfo()## R version 3.1.2 (2014-10-31)
## Platform: i686-pc-linux-gnu (32-bit)
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=pt_BR.UTF-8
## [4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods splines grid stats graphics grDevices utils datasets
## [9] base
##
## other attached packages:
## [1] wzRfun_0.4 ellipse_0.3-8 multcomp_1.3-7 TH.data_1.0-5
## [5] mvtnorm_1.0-1 doBy_4.5-12 survival_2.37-7 reshape_0.8.5
## [9] plyr_1.8.1 alr3_2.0.5 car_2.0-22 gridExtra_0.9.1
## [13] latticeExtra_0.6-26 RColorBrewer_1.0-5 lattice_0.20-29 rmarkdown_0.3.3
## [17] knitr_1.8
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.4 evaluate_0.5.5 formatR_1.0 htmltools_0.2.6 MASS_7.3-35
## [6] Matrix_1.1-4 nnet_7.3-8 Rcpp_0.11.3 sandwich_2.3-2 stringr_0.6.2
## [11] tools_3.1.2 yaml_2.1.13 zoo_1.7-11