source("config/setup.R") detach("package:ggplot2")
#----------------------------------------------------------------------- # 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) # 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 p <- 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.3
O 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))) summary(da) # 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))
# Frequencia de cada tipo de censura. xtabs(~status, data = db) # Áreas abaixo da curva da altura e do diâmetro. xyplot(aacalt + aacdac ~ dos | reorder(gen, aacdac, sum), data = db, xlab = "Dose de cálcio (codificada)", ylab = "Área abaixo da curva", auto.key = TRUE, type = c("p", "a")) # subset(db, status == 2) # subset(da, gen == "8" & dos == "1")
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)) # 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) # 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) # Reduzindo o modelo pela retirada da interação. m0 <- update(m0, . ~ gen + dosep) anova(m0) # Estimativas dos parâmetros. summary(m0)
Nã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 # 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 # Atribui nomes às linhas da matriz. L <- lsm$L rownames(L) <- levels(dc$gen) # Obtém os contrastes dois a dois. grid <- apmc(X = L, model = m0, focus = "gen", test = "fdr", cld2 = FALSE) grid <- grid[order(grid$fit, decreasing = TRUE), ] ocld <- ordered_cld(grid$cld, grid$fit) grid$cld <- ocld # x <- attr(ocld, "ind") * 1 # print.table(local({x[x == 0] <- NA; x})) # 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)", # par.settings = list(layout.widths = list(right.padding = 10)), 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) }) # x <- attr(ocld, "ind") # index <- which(x[nrow(x):1, ], arr.ind = TRUE) # trellis.focus("panel", column = 1, row = 1, clip.off = TRUE) # xcor <- 1.03 + (index[, 2] - 1)/50 # grid.segments(x0 = unit(xcor, "npc"), # x1 = unit(xcor, "npc"), # y0 = unit(index[, 1] + 0.5, units = "native"), # y1 = unit(index[, 1] - 0.5, units = "native"), # gp = gpar(lwd = 2, col = "red")) # trellis.unfocus()
# 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) # 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) # Reduzindo o modelo pela retirada da interação. m0 <- update(m0, . ~ gen + dosep) anova(m0) # Estimativas dos parâmetros. summary(m0)
Nã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 # 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 # Atribui nomes às linhas da matriz. L <- lsm$L rownames(L) <- levels(dc$gen) # Obtém os contrastes dois a dois. grid <- apmc(X = L, model = m0, focus = "gen", test = "fdr", cld2 = FALSE) grid <- grid[order(grid$fit, decreasing = TRUE), ] ocld <- ordered_cld(grid$cld, grid$fit) grid$cld <- ocld # 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) })
# 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) # Reduzindo o modelo mantendo apenas genótipo. m0 <- update(m0, . ~ gen) anova(m0) # Estimativas dos parâmetros. summary(m0)
# Médias ajustadas, considerando um valor fixo de cálcio. lsm <- LSmeans(m0, effect = "gen") lsm # Atribui nomes às linhas da matriz. L <- lsm$L rownames(L) <- levels(dc$gen) # Obtém os contrastes dois a dois. grid <- apmc(X = L, model = m0, focus = "gen", test = "fdr", cld2 = FALSE) grid <- grid[order(grid$fit, decreasing = TRUE), ] ocld <- ordered_cld(grid$cld, grid$fit) grid$cld <- ocld # 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) })
# Utilizar como variável concomitante a última data em foi medida em uma # ancova. str(da) # 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) 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) # Reduzindo o modelo com o abandono do termo quadrático e interação. m0 <- update(m0, . ~ dias + I(dias^2) + gen + dosep) anova(m0) # Estimativas dos parâmetros. summary(m0)
Pelo 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 # 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 # Atribui nomes às linhas da matriz. L <- lsm$L rownames(L) <- levels(dc$gen) # Obtém os contrastes dois a dois. grid <- apmc(X = L, model = m0, focus = "gen", test = "fdr", cld2 = FALSE) grid <- grid[order(grid$fit, decreasing = TRUE), ] ocld <- ordered_cld(grid$cld, grid$fit) grid$cld <- ocld # 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) })
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) # Reduzindo o modelo com o abandono do termo quadrático e interação. m0 <- update(m0, . ~ dias + I(dias^2) + gen) anova(m0) # Estimativas dos parâmetros. summary(m0)
# Ú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 # Atribui nomes às linhas da matriz. L <- lsm$L rownames(L) <- levels(dc$gen) # Obtém os contrastes dois a dois. grid <- apmc(X = L, model = m0, focus = "gen", test = "fdr", cld2 = FALSE) grid <- grid[order(grid$fit, decreasing = TRUE), ] ocld <- ordered_cld(grid$cld, grid$fit) grid$cld <- ocld # 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) })
cat(format(Sys.time(), format = "Atualizado em %d de %B de %Y.\n\n")) sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.