|
Estatística Aplicada à Ciência do Solo
|
Desenvolvimento de mudas de teca em resposta à doses de Cálcio
16 de abril de 2017
Descrição e Análise Exploratória
#-----------------------------------------------------------------------
# Carrega os pacotes necessários.
rm(list = ls())
library(lattice)
library(latticeExtra)
library(grid)
library(gridExtra)
library(plyr)
library(doBy)
library(multcomp)
library(wzRfun)
library(EACS)
# Obtém o intervalo em que a muda morreu.
deathint <- function(data, timevar, respvar) {
o <- order(data[, timevar])
x <- data[o, timevar]
y <- data[o, respvar]
nna <- !is.na(y)
if (sum(nna)) {
i <- which(nna)
imax <- max(i)
st <- 3 # Censura intervalar.
int <- x[imax + 0:1]
} else {
st <- 2 # Censura a esquerda.
int <- rep(min(x), 2)
}
if (is.na(int[2])) {
int[2] <- int[1]
st <- 0 # Censura a direita.
}
r <- c(int, st)
names(r) <- c("left", "right", "status")
return(r)
}Os gráficos abaixo exibem a curva de crescimento das mudas em altura e diâmetro à altura do colo para cada um dos genótipos (paineis) e cada dose de cálcio (cores das linhas). Verifica-se que existem mudas com interrupção das medidas o que é resultado da morte das mudas.
#-----------------------------------------------------------------------
# Análise gráfica exploratória para altura e diâmetro.
# Estrutura dos dados.
str(gen_teca)## 'data.frame': 5800 obs. of 7 variables:
## $ gen : Factor w/ 29 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dose: num 0 0.0184 0.0368 0.0736 0.1472 ...
## $ data: Date, format: "2016-06-15" ...
## $ dac : num 38.6 40.2 43 40 39.8 23.8 33.9 37 40 41.8 ...
## $ alt : num 101 61.7 130 128 128 126 132 121 129 140 ...
## $ mspa: num NA NA NA NA NA NA NA NA NA NA ...
## $ msr : num NA NA NA NA NA NA NA NA NA NA ...# Nomes curtos são mais fáceis de manipular.
da <- gen_teca
# Supondo que a primeira medida foi no dia 2.
da$dias <- as.numeric(da$data - min(da$data) + 2)
xyplot(alt ~ dias | gen,
data = da,
groups = dose,
type = "l",
xlab = "Tempo após a primeira avaliação (dias)",
ylab = "Altura das mudas de teca (mm)",
as.table = TRUE,
scales = list(y = list(log = FALSE)))
xyplot(dac ~ dias | gen,
data = da,
groups = dose,
type = "l",
xlab = "Tempo após a primeira avaliação (dias)",
ylab = "Diâmetro à altura do colo das mudas de teca (mm)",
as.table = TRUE)
Para utilizar o fator métrico dose de cálcio nas análises (dose), é recomendável que se utilize uma transformação na escala das doses de forma a obter um espaçamento mais equidistante entre níveis, o que evita problemas de alancagem para ajustes de modelos de regressão. Em vista disso, definiu-se como critério minimizar a variância entre as diferenças de doses consecutivas na escala unitária por meio de uma função potência das doses transformadas. Na escala unitária, a menor é 0 e a maior é 1. O valor obtido para transformar a dose foi de 0.3.
#-----------------------------------------------------------------------
# Trabalhando com a dose.
# Criando dose categórico.
x <- unique(sort(da$dose))
da$dos <- factor(da$dose,
levels = x,
labels = seq_along(x) - 1)
# Variância das distâncias entre níveis em escala unitária.
esp <- function(p) {
u <- x^p
u <- (u - min(u))
u <- u/max(u)
var(diff(u))
}
# Otimiza para obter o valor de potência para maior uniformidade.
op <- optimize(f = esp, interval = c(0, 1))
op$minimum## [1] 0.3064936p <- seq(0, 1, by = 0.01)
v <- sapply(p, esp)
plot(log(v) ~ p, type = "o")
abline(v = op$minimum)
da$dosep <- da$dose^0.3O gráfico de segmentos a seguir exibe os intervalos de tempo em que as mudas morreram no experimento. A morte das mudas, até onde se sabe, pode não ter relação com os fatores estudados pois depende muito de fatores agronômicos que representam o estado inicial da muda. No entanto, embora o emperimento não tenha sido realizado para estudar fatores relacionados ao “pegamento” das mudas, pode-se verificar se existe influência dos fatores estudados na sobreviênvia das mudas.
#-----------------------------------------------------------------------
# Determinar o intervalo em que a muda morreu.
# Cria a variável unidade experimental.
da$ue <- with(da, interaction(gen, dos, drop = TRUE))
# Número de níveis dos fatores.
with(da, c(gen = nlevels(gen),
dos = nlevels(dos),
ue = nlevels(ue)))## gen dos ue
## 29 10 290summary(da)## gen dose data
## 1 : 200 Min. :0.0000 Min. :2016-06-15
## 2 : 200 1st Qu.:0.0368 1st Qu.:2016-07-18
## 3 : 200 Median :0.2208 Median :2016-08-20
## 4 : 200 Mean :0.9402 Mean :2016-08-20
## 5 : 200 3rd Qu.:1.1776 3rd Qu.:2016-09-22
## 6 : 200 Max. :4.7104 Max. :2016-10-26
## (Other):4600
## dac alt mspa msr
## Min. : 23.30 Min. : 17.0 Min. : 2.96 Min. : 3.07
## 1st Qu.: 52.00 1st Qu.:106.0 1st Qu.:19.36 1st Qu.: 27.81
## Median : 62.00 Median :130.0 Median :36.83 Median : 40.01
## Mean : 67.92 Mean :141.5 Mean :35.76 Mean : 43.79
## 3rd Qu.: 79.03 3rd Qu.:157.0 3rd Qu.:48.67 3rd Qu.: 59.72
## Max. :149.70 Max. :528.0 Max. :93.43 Max. :133.07
## NA's :1728 NA's :1730 NA's :5643 NA's :5643
## dias dos dosep ue
## Min. : 2.00 0 : 580 Min. :0.0000 1.0 : 20
## 1st Qu.: 35.25 1 : 580 1st Qu.:0.3713 2.0 : 20
## Median : 68.50 2 : 580 Median :0.6279 3.0 : 20
## Mean : 68.50 3 : 580 Mean :0.7174 4.0 : 20
## 3rd Qu.:101.75 4 : 580 3rd Qu.:1.0503 5.0 : 20
## Max. :135.00 5 : 580 Max. :1.5919 6.0 : 20
## (Other):2320 (Other):5680# Obtém os intervalos em que as mudas morrem.
db <- ddply(da,
~gen + dose + dos + dosep + ue,
deathint,
timevar = "dias",
respvar = "alt")
# Áreas abaixo da curva (usando semanas e cm como unidades).
dc <- ddply(da,
~gen + dose + dos + dosep + ue,
summarise,
aacalt = auc(dias/7, alt/10),
aacdac = auc(dias/7, dac/10))
# Junta os dois.
db <- merge(db, dc)
# str(db)
# Gráfico de segmentos que indica o intervalo em que a muda morreu.
segplot(dos ~ left + right | gen,
data = db,
as.table = TRUE,
xlab = "Tempo após a primeira avaliação",
ylab = "Dose (codificada)") # +
# layer(panel.text(x = x[subscripts],
# y = z[subscripts],
# labels = db$status[subscripts],
# pos = 2,
# cex = 0.6))Massa seca de parte áerea
A análise dos dados de massa seca de parte aérea foi feita apenas considerando genótipos com 3 ou mais registros, que correspondem ao número de doses de cálcio que restaram ao final do experimento.
dc <- subset(da, data == max(data))
dc <- dc[complete.cases(dc), ]
# Tabela de frequência do número de observações por cela experimental.
addmargins(xtabs(~gen + dos, data = dc))## dos
## gen 0 1 2 3 4 5 6 7 8 9 Sum
## 1 1 1 1 0 0 1 1 1 1 0 7
## 2 1 0 1 1 1 1 1 1 0 0 7
## 3 1 1 1 0 1 1 0 1 1 0 7
## 4 1 0 0 0 0 0 0 1 1 0 3
## 5 0 1 1 0 0 0 0 0 0 0 2
## 6 0 0 0 1 0 0 1 0 0 1 3
## 7 0 0 0 1 0 0 1 1 0 0 3
## 8 0 0 0 0 1 0 0 0 0 0 1
## 9 1 0 0 0 0 1 0 1 0 0 3
## 10 1 0 1 0 1 1 1 1 0 0 6
## 11 0 1 0 1 0 0 1 1 0 0 4
## 12 0 1 1 1 1 1 1 1 0 1 8
## 13 1 1 1 1 0 1 0 1 1 0 7
## 14 1 1 1 0 0 1 1 0 0 0 5
## 15 1 0 0 0 0 1 1 1 1 0 5
## 16 1 1 1 0 0 1 0 1 1 0 6
## 17 1 1 1 0 0 1 1 0 1 0 6
## 18 1 1 0 1 1 0 0 1 0 1 6
## 19 1 1 1 0 1 1 0 0 0 1 6
## 20 1 1 1 0 0 1 0 0 0 1 5
## 21 0 1 1 0 1 1 1 0 1 0 6
## 22 1 1 1 1 1 1 1 0 1 0 8
## 23 1 1 0 0 1 0 1 1 1 0 6
## 24 0 1 1 1 0 1 0 1 0 1 6
## 25 1 1 1 0 1 1 0 1 0 0 6
## 26 1 1 0 1 1 1 0 1 0 0 6
## 27 0 1 1 0 0 1 0 0 1 0 4
## 28 1 1 1 1 1 1 0 1 0 0 7
## 29 1 0 1 0 0 0 0 0 0 1 3
## Sum 20 20 19 11 13 20 13 18 11 7 152# Manter apenas genótipos com 3 ou mais registros.
k <- which(xtabs(~gen, data = dc) >= 3)
# Filtra para os genótipos com 3 ou mais registros.
dc <- droplevels(subset(dc, gen %in% names(k)))
xy1 <- xyplot(mspa + msr ~ dosep | gen,
outer = FALSE,
data = dc,
ylab = "Massa seca (g)",
xlab = "Dose de Ca (escala quase equidistante)",
auto.key = TRUE,
type = "o")
xy2 <- xyplot(mspa + msr ~ dosep | gen,
outer = FALSE,
data = dc,
ylab = "Massa seca (g)",
xlab = "Dose de Ca (escala quase equidistante)",
auto.key = TRUE,
type = c("p", "r"))
# grid.arrange(xy1, xy2, nrow = 1)
xy2
Foi considerado um modelo fatorial duplo contendo os efeitos do fator genótipo (níveis categóricos), do fator dose de cálcio (níveis métricos, expresso na escala transformada) e a interação entre estes fatores. Para representar o efeito de cálcio, foram utilizados polinômio de segundo grau, sendo o grau do polinômio aumentado ou diminuido conforme necessidade de melhorar ou simplificar o ajuste.
# Modelo de efeitos quadráticos por genótipo.
m0 <- lm(mspa ~ gen * (dosep + I(dosep^2)), data = dc)
# Verificação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 22105.4 850.21 2.8684 0.000272 ***
## dosep 1 1477.5 1477.53 4.9848 0.028866 *
## I(dosep^2) 1 27.4 27.37 0.0923 0.762168
## gen:dosep 26 10089.5 388.06 1.3092 0.187858
## gen:I(dosep^2) 26 4301.5 165.44 0.5582 0.950194
## Residuals 68 20155.5 296.41
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo com a remoção dos termos quadráticos.
m0 <- update(m0, . ~ gen * dosep)
# par(mfrow = c(2, 2))
# plot(m0); layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 22105.4 850.21 3.2839 1.278e-05 ***
## dosep 1 1477.5 1477.53 5.7069 0.01887 *
## gen:dosep 26 9978.3 383.78 1.4823 0.08749 .
## Residuals 95 24595.6 258.90
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo pela retirada da interação.
m0 <- update(m0, . ~ gen + dosep)
anova(m0)## Analysis of Variance Table
##
## Response: mspa
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 22105 850.21 2.9755 3.03e-05 ***
## dosep 1 1478 1477.53 5.1710 0.02473 *
## Residuals 121 34574 285.73
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativas dos parâmetros.
summary(m0)##
## Call:
## lm(formula = mspa ~ gen + dosep, data = dc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.957 -12.348 0.087 9.066 48.438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 30.26998 6.75309 4.482 1.69e-05 ***
## gen2 -0.03656 9.03960 -0.004 0.99678
## gen3 25.61210 9.03647 2.834 0.00538 **
## gen4 -9.04181 11.67273 -0.775 0.44008
## gen6 -7.09007 11.71266 -0.605 0.54609
## gen7 3.04863 11.67346 0.261 0.79441
## gen9 19.94146 11.66706 1.709 0.08998 .
## gen10 3.59749 9.40674 0.382 0.70281
## gen11 22.09697 10.59505 2.086 0.03912 *
## gen12 23.75822 8.75296 2.714 0.00761 **
## gen13 28.39832 9.03739 3.142 0.00211 **
## gen14 -7.59291 9.92236 -0.765 0.44562
## gen15 15.22280 9.90683 1.537 0.12700
## gen16 23.38485 9.40502 2.486 0.01427 *
## gen17 18.50569 9.40698 1.967 0.05145 .
## gen18 7.21872 9.40439 0.768 0.44423
## gen19 7.13165 9.40687 0.758 0.44985
## gen20 3.45686 9.89985 0.349 0.72756
## gen21 -9.50005 9.40479 -1.010 0.31445
## gen22 21.23611 8.75317 2.426 0.01674 *
## gen23 -11.64449 9.40472 -1.238 0.21806
## gen24 15.27988 9.40946 1.624 0.10700
## gen25 20.01100 9.41877 2.125 0.03566 *
## gen26 12.43185 9.41623 1.320 0.18924
## gen27 -5.59920 10.59503 -0.528 0.59814
## gen28 8.36645 9.05152 0.924 0.35716
## gen29 9.62461 11.66465 0.825 0.41093
## dosep -7.63221 3.35632 -2.274 0.02473 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.9 on 121 degrees of freedom
## Multiple R-squared: 0.4055, Adjusted R-squared: 0.2729
## F-statistic: 3.057 on 27 and 121 DF, p-value: 1.542e-05Não houve interação entre genótipo e dose de cálcio. Em função disso, o modelo obtido para a predição contém apenas os efeitos aditivos dos fatores. Foi necessário apenas o polinômio de grau um para representar o efeito de cálcio. Os resultados serão representados com as retas ajustadas para cada genótipo.
# Cria um conjunto de valores para predição.
pred <- with(dc,
expand.grid(# dosep = eseq(dose, f = 0)^0.3)
dosep = eseq(dosep, f = 0),
gen = levels(gen)))
ic <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, as.data.frame(ic))
xyplot(mspa ~ dosep | gen,
data = dc,
# type = c("p", "r"),
xlab = "Dose de cálcio (escala transformada)",
ylab = "Massa seca de parte áerea (g)") +
as.layer(xyplot(fit ~ dosep | gen,
data = pred,
type = "l",
ly = pred$lwr,
uy = pred$upr,
cty = "bands",
prepanel = prepanel.cbH,
panel = panel.cbH))
# Médias ajustadas, considerando um valor fixo de cálcio.
lsm <- LSmeans(m0, effect = "gen")
lsm## estimate se df t.stat p.value gen dosep
## 1 25.48490 6.389534 121 3.988538 1.142593e-04 1 0.6269589
## 2 25.44834 6.391887 121 3.981350 1.173558e-04 2 0.6269589
## 3 51.09700 6.389238 121 7.997355 8.598656e-13 3 0.6269589
## 4 16.44309 9.773049 121 1.682493 9.505201e-02 4 0.6269589
## 5 18.39484 9.826002 121 1.872057 6.361290e-02 6 0.6269589
## 6 28.53353 9.774087 121 2.919304 4.184812e-03 7 0.6269589
## 7 45.42637 9.760563 121 4.654072 8.399081e-06 9 0.6269589
## 8 29.08239 6.902114 121 4.213549 4.868040e-05 10 0.6269589
## 9 47.58187 8.452834 121 5.629103 1.194137e-07 11 0.6269589
## 10 49.24313 5.987374 121 8.224495 2.562139e-13 12 0.6269589
## 11 53.88322 6.389882 121 8.432584 8.386814e-14 13 0.6269589
## 12 17.89199 7.584515 121 2.359016 1.992531e-02 14 0.6269589
## 13 40.70770 7.576505 121 5.372887 3.816902e-07 15 0.6269589
## 14 48.86975 6.900963 121 7.081584 1.017314e-10 16 0.6269589
## 15 43.99059 6.902315 121 6.373310 3.501766e-09 17 0.6269589
## 16 32.70363 6.901826 121 4.738402 5.925657e-06 18 0.6269589
## 17 32.61655 6.902219 121 4.725517 6.251622e-06 19 0.6269589
## 18 28.94176 7.560492 121 3.828027 2.059866e-04 20 0.6269589
## 19 15.98486 6.903122 121 2.315598 2.226407e-02 21 0.6269589
## 20 46.72101 5.979800 121 7.813139 2.279566e-12 22 0.6269589
## 21 13.84041 6.902928 121 2.005006 4.719430e-02 23 0.6269589
## 22 40.76478 6.912110 121 5.897589 3.427676e-08 24 0.6269589
## 23 45.49590 6.914780 121 6.579516 1.270363e-09 25 0.6269589
## 24 37.91675 6.911905 121 5.485716 2.296426e-07 26 0.6269589
## 25 19.88570 8.452793 121 2.352560 2.025870e-02 27 0.6269589
## 26 33.85135 6.405298 121 5.284898 5.649573e-07 28 0.6269589
## 27 35.10951 9.759782 121 3.597367 4.666451e-04 29 0.6269589# Dose usada para obter as médias (é a média de cálcio dos dados
# observados).
d <- lsm$grid$dosep[1]^(1/0.3)
d## [1] 0.210926# Atribui nomes às linhas da matriz.
L <- lsm$K
rownames(L) <- levels(dc$gen)
# Obtém os contrastes dois a dois.
grid <- apmc(X = L,
model = m0,
focus = "gen",
test = "fdr",
cld2 = TRUE)
# Gráfico de segmentos com as letras resumo das comparações múltiplas.
segplot(reorder(gen, fit) ~ lwr + upr,
centers = fit,
data = grid,
draw = FALSE,
xlab = "Massa seca de parte áerea (g)",
ylab = "Genótipos (ordenados pela média)",
sub = list(
sprintf("Estimativas correspondentes à %0.2f de Ca", d),
font = 3, cex = 0.8),
txt = sprintf("%0.1f %s", grid$fit, grid$cld)) +
layer({
x <- unit(centers, "native") - unit(0.005, "snpc")
y <- unit(z, "native")
grid.rect(x = x, y = y,
width = Reduce(unit.c,
lapply(txt,
FUN = unit,
x = 1.05,
units = "strwidth")) +
unit(0.005, "snpc"),
height = unit(1, "lines"),
just = "right",
gp = gpar(fill = "white",
col = NA,
fontsize = 12))
grid.text(label = txt,
just = "right",
x = x - unit(0.005, "snpc"),
y = y)
})
Massa seca de raízes
# Modelo de efeitos quadráticos por genótipo.
m0 <- lm(msr ~ gen * (dosep + I(dosep^2)), data = dc)
# Verificação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: msr
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 23001 884.6 1.7046 0.041564 *
## dosep 1 3884 3884.1 7.4840 0.007933 **
## I(dosep^2) 1 377 376.8 0.7260 0.397171
## gen:dosep 26 6946 267.2 0.5148 0.969154
## gen:I(dosep^2) 26 12463 479.3 0.9236 0.576167
## Residuals 68 35291 519.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo com a remoção dos termos quadráticos.
m0 <- update(m0, . ~ gen * dosep)
# par(mfrow = c(2, 2))
# plot(m0); layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: msr
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 23001 884.6 1.7598 0.02569 *
## dosep 1 3884 3884.1 7.7265 0.00656 **
## gen:dosep 26 7320 281.5 0.5601 0.95354
## Residuals 95 47757 502.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo pela retirada da interação.
m0 <- update(m0, . ~ gen + dosep)
anova(m0)## Analysis of Variance Table
##
## Response: msr
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 23001 884.6 1.9435 0.008623 **
## dosep 1 3884 3884.1 8.5331 0.004162 **
## Residuals 121 55077 455.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativas dos parâmetros.
summary(m0)##
## Call:
## lm(formula = msr ~ gen + dosep, data = dc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.807 -11.047 -1.707 9.904 82.182
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 53.0222 8.5234 6.221 7.34e-09 ***
## gen2 -17.9073 11.4093 -1.570 0.11914
## gen3 4.1755 11.4054 0.366 0.71493
## gen4 -18.2065 14.7327 -1.236 0.21893
## gen6 -11.4613 14.7831 -0.775 0.43968
## gen7 4.5578 14.7336 0.309 0.75759
## gen9 3.8148 14.7256 0.259 0.79603
## gen10 -7.8243 11.8727 -0.659 0.51114
## gen11 -12.0316 13.3725 -0.900 0.37005
## gen12 1.9572 11.0475 0.177 0.85967
## gen13 -4.2971 11.4065 -0.377 0.70704
## gen14 15.6134 12.5235 1.247 0.21490
## gen15 10.2654 12.5039 0.821 0.41327
## gen16 6.8593 11.8705 0.578 0.56444
## gen17 7.5859 11.8730 0.639 0.52408
## gen18 -7.5323 11.8697 -0.635 0.52690
## gen19 5.0121 11.8729 0.422 0.67367
## gen20 6.3360 12.4951 0.507 0.61302
## gen21 -12.1416 11.8702 -1.023 0.30841
## gen22 17.4641 11.0478 1.581 0.11654
## gen23 -32.8454 11.8702 -2.767 0.00655 **
## gen24 8.9385 11.8761 0.753 0.45312
## gen25 8.9716 11.8879 0.755 0.45190
## gen26 0.5503 11.8847 0.046 0.96314
## gen27 -22.2116 13.3725 -1.661 0.09930 .
## gen28 -14.7165 11.4244 -1.288 0.20015
## gen29 -6.1375 14.7225 -0.417 0.67751
## dosep -12.3745 4.2362 -2.921 0.00416 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.33 on 121 degrees of freedom
## Multiple R-squared: 0.328, Adjusted R-squared: 0.1781
## F-statistic: 2.188 on 27 and 121 DF, p-value: 0.002119Não houve interação entre genótipo e dose de cálcio. Em função disso, o modelo obtido para a predição contém apenas os efeitos aditivos dos fatores. Foi necessário apenas o polinômio de grau um para representar o efeito de cálcio. Os resultados serão representados com as retas ajustadas para cada genótipo.
# Cria um conjunto de valores para predição.
pred <- with(dc,
expand.grid(# dosep = eseq(dose, f = 0)^0.3)
dosep = eseq(dosep, f = 0),
gen = levels(gen)))
ic <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, as.data.frame(ic))
xyplot(mspa ~ dosep | gen,
data = dc,
# type = c("p", "r"),
xlab = "Dose de cálcio (escala transformada)",
ylab = "Massa seca de raízes (g)") +
as.layer(xyplot(fit ~ dosep | gen,
data = pred,
type = "l",
ly = pred$lwr,
uy = pred$upr,
cty = "bands",
prepanel = prepanel.cbH,
panel = panel.cbH))
# Médias ajustadas, considerando um valor fixo de cálcio.
lsm <- LSmeans(m0, effect = "gen")
lsm## estimate se df t.stat p.value gen dosep
## 1 45.26388 8.064544 121 5.612702 1.287469e-07 1 0.6269589
## 2 27.35658 8.067514 121 3.390955 9.413101e-04 2 0.6269589
## 3 49.43934 8.064170 121 6.130741 1.132526e-08 3 0.6269589
## 4 27.05740 12.335045 121 2.193539 3.017992e-02 4 0.6269589
## 5 33.80255 12.401879 121 2.725599 7.371319e-03 6 0.6269589
## 6 49.82171 12.336355 121 4.038609 9.475646e-05 7 0.6269589
## 7 49.07870 12.319285 121 3.983892 1.162518e-04 9 0.6269589
## 8 37.43955 8.711496 121 4.297717 3.510162e-05 10 0.6269589
## 9 33.23227 10.668736 121 3.114921 2.297682e-03 11 0.6269589
## 10 47.22113 7.556958 121 6.248696 6.416797e-09 12 0.6269589
## 11 40.96683 8.064984 121 5.079592 1.390392e-06 13 0.6269589
## 12 60.87733 9.572789 121 6.359414 3.747442e-09 14 0.6269589
## 13 55.52931 9.562678 121 5.806879 5.243149e-08 15 0.6269589
## 14 52.12317 8.710044 121 5.984260 2.276596e-08 16 0.6269589
## 15 52.84983 8.711750 121 6.066499 1.539849e-08 17 0.6269589
## 16 37.73162 8.711133 121 4.331426 3.075640e-05 18 0.6269589
## 17 50.27594 8.711629 121 5.771129 6.193545e-08 19 0.6269589
## 18 51.59985 9.542468 121 5.407391 3.269640e-07 20 0.6269589
## 19 33.12227 8.712769 121 3.801578 2.266351e-04 21 0.6269589
## 20 62.72801 7.547399 121 8.311209 1.610107e-13 22 0.6269589
## 21 12.41845 8.712524 121 1.425356 1.566289e-01 23 0.6269589
## 22 54.20243 8.724113 121 6.212944 7.626683e-09 24 0.6269589
## 23 54.23551 8.727483 121 6.214336 7.575621e-09 25 0.6269589
## 24 45.81420 8.723855 121 5.251600 6.547130e-07 26 0.6269589
## 25 23.05228 10.668684 121 2.160742 3.268590e-02 27 0.6269589
## 26 30.54741 8.084441 121 3.778543 2.462079e-04 28 0.6269589
## 27 39.12636 12.318300 121 3.176280 1.893023e-03 29 0.6269589# Dose usada para obter as médias (é a média de cálcio dos dados
# observados).
d <- lsm$grid$dosep[1]^(1/0.3)
d## [1] 0.210926# Atribui nomes às linhas da matriz.
L <- lsm$K
rownames(L) <- levels(dc$gen)
# Obtém os contrastes dois a dois.
grid <- apmc(X = L,
model = m0,
focus = "gen",
test = "fdr",
cld2 = TRUE)
# Gráfico de segmentos com as letras resumo das comparações múltiplas.
segplot(reorder(gen, fit) ~ lwr + upr,
centers = fit,
data = grid,
draw = FALSE,
xlab = "Massa seca de raízes (g)",
ylab = "Genótipos (ordenados pela média)",
sub = list(
sprintf("Estimativas correspondentes à %0.2f de Ca", d),
font = 3, cex = 0.8),
txt = sprintf("%0.1f %s", grid$fit, grid$cld)) +
layer({
x <- unit(centers, "native") - unit(0.005, "snpc")
y <- unit(z, "native")
grid.rect(x = x, y = y,
width = Reduce(unit.c,
lapply(txt,
FUN = unit,
x = 1.05,
units = "strwidth")) +
unit(0.005, "snpc"),
height = unit(1, "lines"),
just = "right",
gp = gpar(fill = "white",
col = NA,
fontsize = 12))
grid.text(label = txt,
just = "right",
x = x - unit(0.005, "snpc"),
y = y)
})
Proporção de massa em raízes
# Cria variável que a porcentagem de massa seca de raízes em relação a
# massa seca total da muda.
dc$r <- with(dc, 100 * msr/(msr + mspa))# Modelo de efeitos quadráticos por genótipo.
m0 <- lm(r ~ gen * (dosep + I(dosep^2)), data = dc)
# Verificação dos pressupostos.
par(mfrow = c(2, 2))
plot(m0); layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: r
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 10846.5 417.17 2.3854 0.00225 **
## dosep 1 7.2 7.18 0.0410 0.84007
## I(dosep^2) 1 325.8 325.76 1.8627 0.17682
## gen:dosep 26 5523.2 212.43 1.2147 0.25788
## gen:I(dosep^2) 26 4777.1 183.73 1.0506 0.42085
## Residuals 68 11892.4 174.89
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo mantendo apenas genótipo.
m0 <- update(m0, . ~ gen)
anova(m0)## Analysis of Variance Table
##
## Response: r
## Df Sum Sq Mean Sq F value Pr(>F)
## gen 26 10846 417.17 2.2594 0.001592 **
## Residuals 122 22526 184.64
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativas dos parâmetros.
summary(m0)##
## Call:
## lm(formula = r ~ gen, data = dc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.368 -6.919 0.295 9.255 28.680
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 63.62133 5.13582 12.388 < 2e-16 ***
## gen2 -11.83644 7.26315 -1.630 0.10575
## gen3 -15.17842 7.26315 -2.090 0.03871 *
## gen4 4.27467 9.37668 0.456 0.64928
## gen6 0.48032 9.37668 0.051 0.95923
## gen7 -7.83977 9.37668 -0.836 0.40474
## gen9 -10.97340 9.37668 -1.170 0.24417
## gen10 1.58445 7.55972 0.210 0.83434
## gen11 -20.34053 8.51679 -2.388 0.01846 *
## gen12 -16.13130 7.03251 -2.294 0.02351 *
## gen13 -23.41823 7.26315 -3.224 0.00162 **
## gen14 13.19129 7.95638 1.658 0.09990 .
## gen15 -9.46325 7.95638 -1.189 0.23660
## gen16 -13.08349 7.55972 -1.731 0.08604 .
## gen17 -8.13620 7.55972 -1.076 0.28394
## gen18 -13.76289 7.55972 -1.821 0.07113 .
## gen19 0.08582 7.55972 0.011 0.99096
## gen20 0.62266 7.95638 0.078 0.93775
## gen21 6.83064 7.55972 0.904 0.36801
## gen22 -3.54485 7.03251 -0.504 0.61512
## gen23 -17.66014 7.55972 -2.336 0.02112 *
## gen24 -6.71211 7.55972 -0.888 0.37635
## gen25 -7.21256 7.55972 -0.954 0.34193
## gen26 -8.66591 7.55972 -1.146 0.25390
## gen27 -2.26186 8.51679 -0.266 0.79101
## gen28 -14.78767 7.26315 -2.036 0.04392 *
## gen29 -10.05089 9.37668 -1.072 0.28588
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.59 on 122 degrees of freedom
## Multiple R-squared: 0.325, Adjusted R-squared: 0.1812
## F-statistic: 2.259 on 26 and 122 DF, p-value: 0.001592# Médias ajustadas, considerando um valor fixo de cálcio.
lsm <- LSmeans(m0, effect = "gen")
lsm## estimate se df t.stat p.value gen
## 1 63.62133 5.135820 122 12.387764 2.551870e-23 1
## 2 51.78488 5.135820 122 10.083080 9.268901e-18 2
## 3 48.44291 5.135820 122 9.432361 3.398118e-16 3
## 4 67.89599 7.845095 122 8.654579 2.404057e-14 4
## 5 64.10164 7.845095 122 8.170920 3.267864e-13 6
## 6 55.78156 7.845095 122 7.110374 8.538101e-11 7
## 7 52.64793 7.845095 122 6.710936 6.455023e-10 9
## 8 65.20578 5.547320 122 11.754465 8.506928e-22 10
## 9 43.28080 6.794051 122 6.370396 3.481879e-09 11
## 10 47.49003 4.804120 122 9.885272 2.777272e-17 12
## 11 40.20310 5.135820 122 7.827981 2.029459e-12 13
## 12 76.81261 6.076784 122 12.640339 6.339157e-24 14
## 13 54.15808 6.076784 122 8.912292 5.905162e-15 15
## 14 50.53784 5.547320 122 9.110317 1.997304e-15 16
## 15 55.48512 5.547320 122 10.002150 1.452478e-17 17
## 16 49.85844 5.547320 122 8.987843 3.906904e-15 18
## 17 63.70715 5.547320 122 11.484312 3.815350e-21 19
## 18 64.24398 6.076784 122 10.572036 6.115431e-19 20
## 19 70.45197 5.547320 122 12.700181 4.560180e-24 21
## 20 60.07648 4.804120 122 12.505200 1.334853e-23 22
## 21 45.96119 5.547320 122 8.285296 1.768879e-13 23
## 22 56.90921 5.547320 122 10.258867 3.490743e-18 24
## 23 56.40877 5.547320 122 10.168653 5.762755e-18 25
## 24 54.95541 5.547320 122 9.906661 2.466745e-17 26
## 25 61.35946 6.794051 122 9.031351 3.078918e-15 27
## 26 48.83365 5.135820 122 9.508443 2.233357e-16 28
## 27 53.57043 7.845095 122 6.828526 3.575988e-10 29# Atribui nomes às linhas da matriz.
L <- lsm$K
rownames(L) <- levels(dc$gen)
# Obtém os contrastes dois a dois.
grid <- apmc(X = L,
model = m0,
focus = "gen",
test = "fdr",
cld2 = TRUE)
# Gráfico de segmentos com as letras resumo das comparações múltiplas.
segplot(reorder(gen, fit) ~ lwr + upr,
centers = fit,
data = grid,
draw = FALSE,
xlab = "Percentual de massa de raízes",
ylab = "Genótipos (ordenados pela média)",
txt = sprintf("%0.1f %s", grid$fit, grid$cld)) +
layer({
x <- unit(centers, "native") - unit(0.005, "snpc")
y <- unit(z, "native")
grid.rect(x = x, y = y,
width = Reduce(unit.c,
lapply(txt,
FUN = unit,
x = 1.05,
units = "strwidth")) +
unit(0.005, "snpc"),
height = unit(1, "lines"),
just = "right",
gp = gpar(fill = "white",
col = NA,
fontsize = 12))
grid.text(label = txt,
just = "right",
x = x - unit(0.005, "snpc"),
y = y)
})
Altura
# Utilizar como variável concomitante a última data em foi medida em uma
# ancova.
str(da)## 'data.frame': 5800 obs. of 11 variables:
## $ gen : Factor w/ 29 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dose : num 0 0.0184 0.0368 0.0736 0.1472 ...
## $ data : Date, format: "2016-06-15" ...
## $ dac : num 38.6 40.2 43 40 39.8 23.8 33.9 37 40 41.8 ...
## $ alt : num 101 61.7 130 128 128 126 132 121 129 140 ...
## $ mspa : num NA NA NA NA NA NA NA NA NA NA ...
## $ msr : num NA NA NA NA NA NA NA NA NA NA ...
## $ dias : num 2 2 2 2 2 2 2 2 2 2 ...
## $ dos : Factor w/ 10 levels "0","1","2","3",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ dosep: num 0 0.302 0.371 0.457 0.563 ...
## $ ue : Factor w/ 290 levels "1.0","2.0","3.0",..: 1 30 59 88 117 146 175 204 233 262 ...# Função que extraí o registro completo mais velho de uma unidade
# experimental (última observação completa registrada).
last <- function(data, time, resp) {
tm <- data[, time]
data <- data[order(tm), ]
i <- max(which(!is.na(data[, resp])))
data[i, ]
}
# Reduz para uma tabela com uma linha por UE.
dc <- ddply(da,
~ue,
last, time = "data", resp = "alt")
str(dc)## 'data.frame': 290 obs. of 11 variables:
## $ gen : Factor w/ 29 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ dose : num 0 0 0 0 0 0 0 0 0 0 ...
## $ data : Date, format: "2016-10-26" ...
## $ dac : num 88 76.1 105.9 51.6 109.2 ...
## $ alt : num 325 182 244 112 400 179 334 60 468 380 ...
## $ mspa : num 42.59 44.47 39.64 4.22 NA ...
## $ msr : num 47.9 53.5 49 11.4 78.2 ...
## $ dias : num 135 135 135 135 135 58 128 30 135 135 ...
## $ dos : Factor w/ 10 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ dosep: num 0 0 0 0 0 0 0 0 0 0 ...
## $ ue : Factor w/ 290 levels "1.0","2.0","3.0",..: 1 2 3 4 5 6 7 8 9 10 ...xyplot(alt + dac ~ dias,
outer = TRUE,
data = da,
xlab = "Tempo após a primeira avaliação (dias)",
ylab = "Altura das mudas de teca (mm)",
col = "gray",
type = c("p", "smooth")) +
as.layer(xyplot(alt + dac ~ I(dias + 1),
data = dc,
outer = TRUE,
col = "blue",
pch = 19,
type = c("p", "smooth")))
Para análise das variáveis altura e diâmetro à altura do colo será considerado análise de covariância. A covariável usada para corrigir os valores observados das respostas será o número de dias no último registro das variáveis. Dessa forma, todas as observações serão aproveitadas na análise, mesmo as que morreram antes do final do experimento.
Verfica-se que uma relação não linear entre a altura (diâmetro) em função do dia em que é medida. Cada ponto em azul no gráfico representa o último valor registrado das variáveis para uma unidade experimental, ou seja, é o valor das variáveis mais próximo do final do experimento, imediatamente antes da muda morrer. Dada a relação quadrática exibida, o efeito do tempo será considerado no modelo com um termo polinomial de grau 2.
# Modelo de efeitos quadráticos por genótipo.
m0 <- lm(alt ~ dias + I(dias^2) + gen * (dosep + I(dosep^2)),
data = dc)
# # Verificação dos pressupostos.
# par(mfrow = c(2, 2))
# plot(m0); layout(1)
# MASS::boxcox(m0)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: alt
## Df Sum Sq Mean Sq F value Pr(>F)
## dias 1 1346945 1346945 310.0027 < 2.2e-16 ***
## I(dias^2) 1 26272 26272 6.0466 0.0147833 *
## gen 28 301451 10766 2.4778 0.0001522 ***
## dosep 1 107996 107996 24.8555 1.336e-06 ***
## I(dosep^2) 1 10278 10278 2.3656 0.1256161
## gen:dosep 28 114378 4085 0.9402 0.5564251
## gen:I(dosep^2) 28 73378 2621 0.6031 0.9435461
## Residuals 200 868989 4345
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo com o abandono do termo quadrático e interação.
m0 <- update(m0, . ~ dias + I(dias^2) + gen + dosep)
anova(m0)## Analysis of Variance Table
##
## Response: alt
## Df Sum Sq Mean Sq F value Pr(>F)
## dias 1 1346945 1346945 324.4210 < 2.2e-16 ***
## I(dias^2) 1 26272 26272 6.3278 0.0125 *
## gen 28 301451 10766 2.5931 4.704e-05 ***
## dosep 1 107996 107996 26.0116 6.595e-07 ***
## Residuals 257 1067024 4152
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativas dos parâmetros.
summary(m0)##
## Call:
## lm(formula = alt ~ dias + I(dias^2) + gen + dosep, data = dc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -152.79 -37.90 -4.71 29.92 214.71
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 145.160655 31.157756 4.659 5.11e-06 ***
## dias -1.026481 0.839520 -1.223 0.22256
## I(dias^2) 0.013895 0.005057 2.747 0.00643 **
## gen2 -24.743539 28.834653 -0.858 0.39163
## gen3 20.672477 28.844344 0.717 0.47422
## gen4 4.959274 29.131254 0.170 0.86496
## gen5 9.890866 29.010954 0.341 0.73343
## gen6 28.156789 29.152219 0.966 0.33503
## gen7 22.478053 29.034382 0.774 0.43953
## gen8 0.041108 30.097746 0.001 0.99891
## gen9 64.167965 28.941270 2.217 0.02749 *
## gen10 12.765646 28.880980 0.442 0.65885
## gen11 14.896021 28.942090 0.515 0.60722
## gen12 27.494495 28.817248 0.954 0.34093
## gen13 75.100000 28.816118 2.606 0.00969 **
## gen14 77.156956 28.959852 2.664 0.00820 **
## gen15 63.562499 28.829519 2.205 0.02836 *
## gen16 44.974263 28.870276 1.558 0.12051
## gen17 19.892636 28.833911 0.690 0.49088
## gen18 21.056149 29.017472 0.726 0.46872
## gen19 40.135602 28.863130 1.391 0.16557
## gen20 32.522052 28.920526 1.125 0.26184
## gen21 14.874796 28.867703 0.515 0.60680
## gen22 34.150331 28.837773 1.184 0.23742
## gen23 -61.810695 29.056100 -2.127 0.03435 *
## gen24 10.487665 28.838070 0.364 0.71640
## gen25 79.385805 28.866006 2.750 0.00638 **
## gen26 83.053311 28.971415 2.867 0.00449 **
## gen27 14.482783 29.224624 0.496 0.62062
## gen28 -9.496780 28.944889 -0.328 0.74310
## gen29 -12.135498 29.147721 -0.416 0.67751
## dosep -42.705333 8.373344 -5.100 6.59e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 64.43 on 257 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.6256, Adjusted R-squared: 0.5804
## F-statistic: 13.85 on 31 and 257 DF, p-value: < 2.2e-16Pelo quadro de ANOVA não existe interação entre genótipo e dose de cálcio. Verifica-se a existência apenas dos efeitos aditivos. No gráfico a seguir é representada a relação entre altura final e dose de cálcio para cada genótipo. Segundo o gráfico e também a estimativa do coeficiente do efeito de Ca, a altura do genótipo diminui com o aumento da dose de cálcio.
# Último dia de medição das mudas.
dmax <- max(dc$dias, na.rm = TRUE)
# Cria um conjunto de valores para predição.
pred <- with(dc,
expand.grid(
dosep = eseq(dosep, f = 0),
dias = dmax,
gen = levels(gen)))
ic <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, as.data.frame(ic))
xyplot(alt ~ dosep | gen,
data = subset(dc, dias == max(dias, na.rm = TRUE)),
# type = c("p", "r"),
xlab = "Dose de cálcio (escala transformada)",
ylab = "Altura das mudas (mm)") +
as.layer(xyplot(fit ~ dosep | gen,
data = pred,
type = "l",
ly = pred$lwr,
uy = pred$upr,
cty = "bands",
prepanel = prepanel.cbH,
panel = panel.cbH))
As médias de altura dos genótipos serão comparadas considerando a altura final, ou seja, aquela prevista para os 135 dias.
# Médias ajustadas, considerando um valor fixo de cálcio para o último
# dia de medição.
lsm <- LSmeans(m0,
effect = "gen",
at = list("dias" = dmax, "I(dias^2)" = dmax^2))
lsm## estimate se df t.stat p.value gen dias I(dias^2)
## 1 229.1299 20.45033 257 11.204214 5.474996e-24 1 135 18225
## 2 204.3863 20.55556 257 9.943114 6.405929e-20 2 135 18225
## 3 249.8023 20.46428 257 12.206749 2.463392e-27 3 135 18225
## 4 234.0891 21.22625 257 11.028286 2.072367e-23 4 135 18225
## 5 239.0207 20.97922 257 11.393213 1.300236e-24 5 135 18225
## 6 257.2867 21.24493 257 12.110495 5.203910e-27 6 135 18225
## 7 251.6079 21.03579 257 11.960946 1.658161e-26 7 135 18225
## 8 229.1710 22.62822 257 10.127665 1.671569e-20 8 135 18225
## 9 293.2978 20.81983 257 14.087427 8.727442e-34 9 135 18225
## 10 241.8955 20.57914 257 11.754403 8.164078e-26 10 135 18225
## 11 244.0259 20.83084 257 11.714643 1.108617e-25 11 135 18225
## 12 256.6244 20.44675 257 12.550861 1.679643e-28 12 135 18225
## 13 304.2299 20.45033 257 14.876526 1.568915e-36 13 135 18225
## 14 306.2868 20.83551 257 14.700230 6.457492e-36 14 135 18225
## 15 292.6924 20.53780 257 14.251400 2.353023e-34 15 135 18225
## 16 274.1041 20.67383 257 13.258510 6.379765e-31 16 135 18225
## 17 249.0225 20.55448 257 12.115244 5.015574e-27 17 135 18225
## 18 250.1860 20.91005 257 11.964871 1.608562e-26 18 135 18225
## 19 269.2655 20.65352 257 13.037266 3.667506e-30 19 135 18225
## 20 261.6519 20.76175 257 12.602597 1.119976e-28 20 135 18225
## 21 244.0047 20.65461 257 11.813567 5.175470e-26 21 135 18225
## 22 263.2802 20.46009 257 12.867987 1.392729e-29 22 135 18225
## 23 167.3192 20.96875 257 7.979454 4.883165e-14 23 135 18225
## 24 239.6175 20.53947 257 11.666201 1.608763e-25 24 135 18225
## 25 308.5157 20.66247 257 14.931208 1.011369e-36 25 135 18225
## 26 312.1832 20.76345 257 15.035225 4.385870e-37 26 135 18225
## 27 243.6126 21.15186 257 11.517314 5.038111e-25 27 135 18225
## 28 219.6331 20.62396 257 10.649414 3.552397e-22 28 135 18225
## 29 216.9944 21.21943 257 10.226211 8.123827e-21 29 135 18225
## dosep
## 1 0.71885
## 2 0.71885
## 3 0.71885
## 4 0.71885
## 5 0.71885
## 6 0.71885
## 7 0.71885
## 8 0.71885
## 9 0.71885
## 10 0.71885
## 11 0.71885
## 12 0.71885
## 13 0.71885
## 14 0.71885
## 15 0.71885
## 16 0.71885
## 17 0.71885
## 18 0.71885
## 19 0.71885
## 20 0.71885
## 21 0.71885
## 22 0.71885
## 23 0.71885
## 24 0.71885
## 25 0.71885
## 26 0.71885
## 27 0.71885
## 28 0.71885
## 29 0.71885# Dose usada para obter as médias (é a média de cálcio dos dados
# observados).
d <- lsm$grid$dosep[1]^(1/0.3)
d## [1] 0.3327573# Atribui nomes às linhas da matriz.
L <- lsm$K
rownames(L) <- levels(dc$gen)
# Obtém os contrastes dois a dois.
grid <- apmc(X = L,
model = m0,
focus = "gen",
test = "fdr",
cld2 = FALSE)
# Gráfico de segmentos com as letras resumo das comparações múltiplas.
segplot(reorder(gen, fit) ~ lwr + upr,
centers = fit,
data = grid,
draw = FALSE,
xlab = sprintf("Altura das mudas aos %d dias (mm)", dmax),
ylab = "Genótipos (ordenados pela média)",
sub = list(
sprintf("Estimativas correspondentes à %0.2f de Ca", d),
font = 3, cex = 0.8),
txt = sprintf("%0.1f %s", grid$fit, grid$cld)) +
layer({
x <- unit(centers, "native") - unit(0.005, "snpc")
y <- unit(z, "native")
grid.rect(x = x, y = y,
width = Reduce(unit.c,
lapply(txt,
FUN = unit,
x = 1.05,
units = "strwidth")) +
unit(0.005, "snpc"),
height = unit(1, "lines"),
just = "right",
gp = gpar(fill = "white",
col = NA,
fontsize = 12))
grid.text(label = txt,
just = "right",
x = x - unit(0.005, "snpc"),
y = y)
})
Diâmetro à altura do colo
A análise para diâmetro à altura do colo segue o mesmo roteiro da análise para altura das mudas.
# Modelo de efeitos quadráticos por genótipo.
m0 <- lm(dac ~ dias + I(dias^2) + gen * (dosep + I(dosep^2)),
data = dc)
# # Verificação dos pressupostos.
# par(mfrow = c(2, 2))
# plot(m0); layout(1)
# MASS::boxcox(m0)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: dac
## Df Sum Sq Mean Sq F value Pr(>F)
## dias 1 140642 140642 512.5683 < 2.2e-16 ***
## I(dias^2) 1 1680 1680 6.1240 0.01417 *
## gen 28 27876 996 3.6284 5.305e-08 ***
## dosep 1 466 466 1.6977 0.19409
## I(dosep^2) 1 231 231 0.8426 0.35977
## gen:dosep 28 10831 387 1.4097 0.09281 .
## gen:I(dosep^2) 28 5488 196 0.7143 0.85411
## Residuals 200 54877 274
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Reduzindo o modelo com o abandono do termo quadrático e interação.
m0 <- update(m0, . ~ dias + I(dias^2) + gen)
anova(m0)## Analysis of Variance Table
##
## Response: dac
## Df Sum Sq Mean Sq F value Pr(>F)
## dias 1 140642 140642 504.7161 < 2.2e-16 ***
## I(dias^2) 1 1680 1680 6.0302 0.01472 *
## gen 28 27876 996 3.5728 3.223e-08 ***
## Residuals 258 71893 279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Estimativas dos parâmetros.
summary(m0)##
## Call:
## lm(formula = dac ~ dias + I(dias^2) + gen, data = dc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -55.793 -9.264 0.755 9.833 49.427
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 56.600728 7.891835 7.172 7.8e-12 ***
## dias -0.263033 0.217329 -1.210 0.22727
## I(dias^2) 0.004164 0.001308 3.183 0.00164 **
## gen2 -8.018830 7.470035 -1.073 0.28406
## gen3 15.496212 7.472634 2.074 0.03910 *
## gen4 -5.283248 7.544041 -0.700 0.48436
## gen5 6.101576 7.514009 0.812 0.41752
## gen6 14.493647 7.549360 1.920 0.05598 .
## gen7 20.382237 7.519841 2.710 0.00717 **
## gen8 4.675344 7.795536 0.600 0.54920
## gen9 12.473500 7.496640 1.664 0.09735 .
## gen10 6.614733 7.481964 0.884 0.37747
## gen11 17.733836 7.496809 2.366 0.01874 *
## gen12 19.657242 7.465616 2.633 0.00897 **
## gen13 21.320000 7.465323 2.856 0.00464 **
## gen14 8.587985 7.501422 1.145 0.25333
## gen15 18.410475 7.468671 2.465 0.01435 *
## gen16 14.878431 7.478959 1.989 0.04772 *
## gen17 14.821973 7.469839 1.984 0.04829 *
## gen18 -5.479486 7.517263 -0.729 0.46671
## gen19 11.124523 7.477072 1.488 0.13802
## gen20 6.664457 7.491505 0.890 0.37451
## gen21 7.794979 7.478214 1.042 0.29822
## gen22 21.090577 7.470931 2.823 0.00513 **
## gen23 -18.818336 7.527289 -2.500 0.01304 *
## gen24 20.405763 7.470871 2.731 0.00674 **
## gen25 17.930494 7.477796 2.398 0.01720 *
## gen26 12.710496 7.504954 1.694 0.09155 .
## gen27 -2.040574 7.569789 -0.270 0.78771
## gen28 8.197417 7.498590 1.093 0.27533
## gen29 3.192171 7.548354 0.423 0.67272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.69 on 258 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.703, Adjusted R-squared: 0.6685
## F-statistic: 20.36 on 30 and 258 DF, p-value: < 2.2e-16# Último dia de medição das mudas.
dmax <- max(dc$dias, na.rm = TRUE)
# Cria um conjunto de valores para predição.
pred <- with(dc,
expand.grid(
dosep = eseq(dosep, f = 0),
dias = dmax,
gen = levels(gen)))
ic <- predict(m0, newdata = pred, interval = "confidence")
pred <- cbind(pred, as.data.frame(ic))
xyplot(dac ~ dosep | gen,
data = subset(dc, dias == max(dias, na.rm = TRUE)),
# type = c("p", "r"),
xlab = "Dose de cálcio (escala transformada)",
ylab = "Altura das mudas (mm)") +
as.layer(xyplot(fit ~ dosep | gen,
data = pred,
type = "l",
ly = pred$lwr,
uy = pred$upr,
cty = "bands",
prepanel = prepanel.cbH,
panel = panel.cbH))
# Médias ajustadas, considerando um valor fixo de cálcio para o último
# dia de medição.
lsm <- LSmeans(m0,
effect = "gen",
at = list("dias" = dmax, "I(dias^2)" = dmax^2))
lsm## estimate se df t.stat p.value gen dias I(dias^2)
## 1 96.97325 5.297346 258 18.30601 1.523758e-48 1 135 18225
## 2 88.95442 5.323903 258 16.70850 5.628797e-43 2 135 18225
## 3 112.46946 5.300884 258 21.21711 1.677498e-58 3 135 18225
## 4 91.69000 5.491145 258 16.69779 6.135638e-43 4 135 18225
## 5 103.07483 5.429358 258 18.98472 6.841776e-51 5 135 18225
## 6 111.46690 5.495773 258 20.28230 2.444830e-55 6 135 18225
## 7 117.35549 5.443533 258 21.55870 1.196763e-59 7 135 18225
## 8 101.64860 5.856792 258 17.35568 3.086408e-45 8 135 18225
## 9 109.44675 5.389529 258 20.30729 2.009991e-55 9 135 18225
## 10 103.58798 5.329710 258 19.43595 1.914550e-52 10 135 18225
## 11 114.70709 5.392285 258 21.27244 1.092894e-58 11 135 18225
## 12 116.63049 5.296405 258 22.02069 3.432179e-61 12 135 18225
## 13 118.29325 5.297346 258 22.33066 3.208126e-62 13 135 18225
## 14 105.56124 5.393521 258 19.57186 6.540590e-53 14 135 18225
## 15 115.38373 5.319152 258 21.69213 4.280895e-60 15 135 18225
## 16 111.85168 5.353490 258 20.89323 2.073085e-57 16 135 18225
## 17 111.79522 5.323604 258 20.99991 9.045598e-58 17 135 18225
## 18 91.49377 5.415242 258 16.89560 1.247987e-43 18 135 18225
## 19 108.09777 5.348126 258 20.21227 4.233173e-55 19 135 18225
## 20 103.63771 5.375036 258 19.28130 6.510263e-52 20 135 18225
## 21 104.76823 5.348289 258 19.58911 5.707617e-53 21 135 18225
## 22 118.06383 5.299790 258 22.27708 4.829067e-62 22 135 18225
## 23 78.15492 5.430518 258 14.39180 7.110685e-35 23 135 18225
## 24 117.37902 5.319521 258 22.06571 2.431064e-61 24 135 18225
## 25 114.90375 5.350385 258 21.47579 2.269046e-59 25 135 18225
## 26 109.68375 5.376145 258 20.40193 9.580199e-56 26 135 18225
## 27 94.93268 5.475013 258 17.33926 3.521553e-45 27 135 18225
## 28 105.17067 5.341611 258 19.68894 2.596020e-53 28 135 18225
## 29 100.16542 5.489465 258 18.24685 2.444286e-48 29 135 18225# Atribui nomes às linhas da matriz.
L <- lsm$K
rownames(L) <- levels(dc$gen)
# Obtém os contrastes dois a dois.
grid <- apmc(X = L,
model = m0,
focus = "gen",
test = "fdr",
cld2 = TRUE)
# Gráfico de segmentos com as letras resumo das comparações múltiplas.
segplot(reorder(gen, fit) ~ lwr + upr,
centers = fit,
data = grid,
draw = FALSE,
xlab = sprintf("Diâmetro à altura do colo aos %d dias (mm)", dmax),
ylab = "Genótipos (ordenados pela média)",
# sub = list(
# sprintf("Estimativas correspondentes à %0.2f de Ca", d),
# font = 3, cex = 0.8),
txt = sprintf("%0.1f %s", grid$fit, grid$cld)) +
layer({
x <- unit(centers, "native") - unit(0.005, "snpc")
y <- unit(z, "native")
grid.rect(x = x, y = y,
width = Reduce(unit.c,
lapply(txt,
FUN = unit,
x = 1.05,
units = "strwidth")) +
unit(0.005, "snpc"),
height = unit(1, "lines"),
just = "right",
gp = gpar(fill = "white",
col = NA,
fontsize = 12))
grid.text(label = txt,
just = "right",
x = x - unit(0.005, "snpc"),
y = y)
})
Informações da Sessão
## Atualizado em 16 de abril de 2017.
##
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=pt_BR.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=pt_BR.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=pt_BR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=pt_BR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] methods grid stats graphics grDevices utils
## [7] datasets base
##
## other attached packages:
## [1] multcomp_1.4-6 TH.data_1.0-8 MASS_7.3-45
## [4] survival_2.40-1 mvtnorm_1.0-5 doBy_4.5-15
## [7] plyr_1.8.4 gridExtra_2.2.1 EACS_0.0-3
## [10] wzRfun_0.77 captioner_2.2.3 latticeExtra_0.6-28
## [13] RColorBrewer_1.1-2 lattice_0.20-34 rmarkdown_1.3
## [16] knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.9 tools_3.3.3 digest_0.6.12
## [4] evaluate_0.10 memoise_1.0.0 gtable_0.2.0
## [7] Matrix_1.2-8 commonmark_1.1 yaml_2.1.14
## [10] rpanel_1.1-3 withr_1.0.2 stringr_1.1.0
## [13] roxygen2_6.0.1 xml2_1.1.1 devtools_1.12.0
## [16] rprojroot_1.2 R6_2.2.0 magrittr_1.5
## [19] backports_1.0.5 codetools_0.2-15 htmltools_0.3.5
## [22] splines_3.3.3 sandwich_2.3-4 stringi_1.1.2
## [25] zoo_1.7-14
Estatística Aplicada à Ciência do Solo |
github.com/walmes/EACS |