# https://github.com/walmes/wzRfun # devtools::install_github("walmes/wzRfun") library(wzRfun) pks <- c("lattice", "gridExtra", "lme4", "lmerTest", "doBy", "multcomp") sapply(pks, library, character.only = TRUE, logical.return = TRUE)
library(RDASC)
source("config/setup.R")
# Object structure. str(sugarcane_straw) # Frequencies. ftable(xtabs(~palha + N + K, data = sugarcane_straw)) # Checking if is a complete cases dataset. all(complete.cases(sugarcane_straw)) # Descriptive measures. summary(sugarcane_straw) # A more detailed description. # Hmisc::describe(sugarcane_straw) # Create factors. sugarcane_straw <- within(sugarcane_straw, { # Convert integer to factor. palha <- factor(sugarcane_straw$palha) nit <- factor(sugarcane_straw$N) pot <- factor(sugarcane_straw$K) # Create a factor to represent main plots. ue <- interaction(bloc, palha) })
#----------------------------------------------------------------------- # To minimize modifications of code when analysing diferent responses. # Define the legends for factors and variables. leg <- list(N = expression("Nitrogênio"~(kg~ha^{-1})), K = expression("Potássio"~(kg~ha^{-1})), y = expression("Número de colmos por metro"), palha = c("Cobertura do solo", "Com palha", "Sem palha")) # Use name y for the response. sugarcane_straw$y <- sugarcane_straw$ncm #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) # plot(pk) plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0) # Drop non relevant terms. # m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue)) m1 <- lmer(y ~ palha * nit + pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m1, m0) anova(m1) # summary(m1) # Reducing by using linear effect for N and K. m2 <- lmer(y ~ palha * N + K + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m2, m0) anova(m2) summary(m2) # Final model. mod <- m2 #----------------------------------------------------------------------- # Fitted values. # Experimental values to get estimates. grid <- with(sugarcane_straw, expand.grid(palha = levels(palha), N = seq(min(N), max(N), length.out = 10), K = seq(min(N), max(N), length.out = 3), KEEP.OUT.ATTRS = FALSE)) # Matrix of fixed effects. X <- model.matrix(terms(mod), data = cbind(grid, y = 0)) # Confidence intervals. ci <- confint(glht(mod, linfct = X), calpha = univariate_calpha())$confint colnames(ci)[1] <- "fit" grid <- cbind(grid, ci) str(grid) # Sample averages. # averages <- aggregate(nobars(formula(mod)), # data = sugarcane_straw, FUN = mean) # # xyplot(y ~ N | factor(K), groups = palha, data = averages, # type = "o", as.table = TRUE) xyplot(fit ~ N | factor(K), groups = palha, data = grid, type = "l", as.table = TRUE, layout = c(NA, 1), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$palha[1], cex.title = 1.1, text = c(leg$palha[-1]), columns = 2, points = FALSE, lines = TRUE), uy = grid$upr, ly = grid$lwr, cty = "bands", alpha = 0.25, fill = "gray50", prepanel = prepanel.cbH, panel.groups = panel.cbH, panel = panel.superpose)
#----------------------------------------------------------------------- leg$y <- "Peso médio de colmo (kg)" sugarcane_straw$y <- sugarcane_straw$pmc #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) # plot(pk) plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0) # Drop non relevant terms. # m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue)) m1 <- lmer(y ~ palha + nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m1, m0) anova(m1) # summary(m1) # Reducing by using linear effect for N and K. m2 <- lmer(y ~ palha + N * K + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m2, m0) anova(m2) summary(m2) # Final model. mod <- m2 #----------------------------------------------------------------------- # Fitted values. # Experimental values to get estimates. grid <- with(sugarcane_straw, expand.grid(palha = levels(palha), N = seq(min(N), max(N), length.out = 10), K = seq(min(N), max(N), length.out = 3), KEEP.OUT.ATTRS = FALSE)) # Matrix of fixed effects. X <- model.matrix(terms(mod), data = cbind(grid, y = 0)) # Confidence intervals. ci <- confint(glht(mod, linfct = X), calpha = univariate_calpha())$confint colnames(ci)[1] <- "fit" grid <- cbind(grid, ci) str(grid) # Sample averages. # averages <- aggregate(nobars(formula(mod)), # data = sugarcane_straw, FUN = mean) # # xyplot(y ~ N | factor(K), groups = palha, data = averages, # type = "o", as.table = TRUE) xyplot(fit ~ N | factor(K), groups = palha, data = grid, type = "l", as.table = TRUE, layout = c(NA, 1), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$palha[1], cex.title = 1.1, text = c(leg$palha[-1]), columns = 2, points = FALSE, lines = TRUE), uy = grid$upr, ly = grid$lwr, cty = "bands", alpha = 0.25, fill = "gray50", prepanel = prepanel.cbH, panel.groups = panel.cbH, panel = panel.superpose)
#----------------------------------------------------------------------- leg$y <- expression("Produção de cana"~(ton~ha^{-1})) sugarcane_straw$y <- sugarcane_straw$tch #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) # plot(pk) plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0) # Drop non relevant terms. # m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue)) m1 <- lmer(y ~ palha + nit + pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m1, m0) anova(m1) # summary(m1) # Reducing by using linear effect for N and K. m2 <- lmer(y ~ palha + N + K + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m2, m0) anova(m2) summary(m2) # Final model. mod <- m2 #----------------------------------------------------------------------- # Fitted values. # Experimental values to get estimates. grid <- with(sugarcane_straw, expand.grid(palha = levels(palha), N = seq(min(N), max(N), length.out = 10), K = seq(min(N), max(N), length.out = 3), KEEP.OUT.ATTRS = FALSE)) # Matrix of fixed effects. X <- model.matrix(terms(mod), data = cbind(grid, y = 0)) # Confidence intervals. ci <- confint(glht(mod, linfct = X), calpha = univariate_calpha())$confint colnames(ci)[1] <- "fit" grid <- cbind(grid, ci) str(grid) # Sample averages. # averages <- aggregate(nobars(formula(mod)), # data = sugarcane_straw, FUN = mean) # # xyplot(y ~ N | factor(K), groups = palha, data = averages, # type = "o", as.table = TRUE) xyplot(fit ~ N | factor(K), groups = palha, data = grid, type = "l", as.table = TRUE, layout = c(NA, 1), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$palha[1], cex.title = 1.1, text = c(leg$palha[-1]), columns = 2, points = FALSE, lines = TRUE), uy = grid$upr, ly = grid$lwr, cty = "bands", alpha = 0.25, fill = "gray50", prepanel = prepanel.cbH, panel.groups = panel.cbH, panel = panel.superpose)
#----------------------------------------------------------------------- leg$y <- "Teor de sacarose aparente" sugarcane_straw$y <- sugarcane_straw$pol #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) # plot(pk) plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0)
#----------------------------------------------------------------------- leg$y <- expression("Sacarose"~(ton~ha^{-1})) sugarcane_straw$y <- sugarcane_straw$tsh #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) # plot(pk) plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0) # Drop non relevant terms. # m1 <- update(m0, formula = . ~ palha * nit + pot + (1 | bloc/ue)) m1 <- lmer(y ~ palha + nit + pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m1, m0) anova(m1) # summary(m1) # Reducing by using linear effect for N and K. m2 <- lmer(y ~ palha + N + K + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Test the reduced model. anova(m2, m0) anova(m2) summary(m2) # Final model. mod <- m2 #----------------------------------------------------------------------- # Fitted values. # Experimental values to get estimates. grid <- with(sugarcane_straw, expand.grid(palha = levels(palha), N = seq(min(N), max(N), length.out = 10), K = seq(min(N), max(N), length.out = 3), KEEP.OUT.ATTRS = FALSE)) # Matrix of fixed effects. X <- model.matrix(terms(mod), data = cbind(grid, y = 0)) # Confidence intervals. ci <- confint(glht(mod, linfct = X), calpha = univariate_calpha())$confint colnames(ci)[1] <- "fit" grid <- cbind(grid, ci) str(grid) # Sample averages. # averages <- aggregate(nobars(formula(mod)), # data = sugarcane_straw, FUN = mean) # # xyplot(y ~ N | factor(K), groups = palha, data = averages, # type = "o", as.table = TRUE) xyplot(fit ~ N | factor(K), groups = palha, data = grid, type = "l", as.table = TRUE, layout = c(NA, 1), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$palha[1], cex.title = 1.1, text = c(leg$palha[-1]), columns = 2, points = FALSE, lines = TRUE), uy = grid$upr, ly = grid$lwr, cty = "bands", alpha = 0.25, fill = "gray50", prepanel = prepanel.cbH, panel.groups = panel.cbH, panel = panel.superpose)
#----------------------------------------------------------------------- leg$y <- "Teor de nitrogênio nas folhas" sugarcane_straw$y <- sugarcane_straw$tfn #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) # plot(pk) plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0) # Final model. mod <- m0 #----------------------------------------------------------------------- # Fitted values. # Experimental values to get estimates. grid <- with(sugarcane_straw, expand.grid(palha = levels(palha), # N = seq(min(N), max(N), length.out = 10), # K = seq(min(N), max(N), length.out = 3), nit = levels(nit), pot = levels(pot), KEEP.OUT.ATTRS = FALSE)) # Matrix of fixed effects. X <- model.matrix(terms(mod), data = cbind(grid, y = 0)) # Confidence intervals. ci <- confint(glht(mod, linfct = X), calpha = univariate_calpha())$confint colnames(ci)[1] <- "fit" grid <- cbind(grid, ci) str(grid) # Sample averages. # averages <- aggregate(nobars(formula(mod)), # data = sugarcane_straw, FUN = mean) # # xyplot(y ~ N | factor(K), groups = palha, data = averages, # type = "o", as.table = TRUE) xyplot(fit ~ pot | nit, groups = palha, data = grid, as.table = TRUE, layout = c(NA, 1), type = "o", xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$palha[1], cex.title = 1.1, text = c(leg$palha[-1]), columns = 2, points = FALSE, lines = TRUE), uy = grid$upr, ly = grid$lwr, desloc = 0.2 * scale(as.integer(grid$palha), scale = FALSE), cty = "bars", length = 0.02, alpha = 0.25, fill = "gray50", prepanel = prepanel.cbH, panel.groups = panel.cbH, panel = panel.superpose)
High order interactions present but no smooth effect of numeric factors.
#----------------------------------------------------------------------- leg$y <- "Teor de potássio nas folhas" sugarcane_straw$y <- sugarcane_straw$tfk #----------------------------------------------------------------------- # Scatter plots. pk <- xyplot(y ~ K | palha, groups = N, data = sugarcane_straw, type = c("p", "a"), xlab = leg$K, ylab = leg$y, auto.key = list(title = leg$N, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) pn <- xyplot(y ~ N | palha, groups = K, data = sugarcane_straw, type = c("p", "a"), xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$K, cex.title = 1.1, columns = 5), strip = strip.custom( strip.names = FALSE, strip.levels = TRUE, factor.levels = c("Com palha", "Sem palha"))) # grid.arrange(pn, pk) plot(pk) # plot(pn)
#----------------------------------------------------------------------- # Model fitting. # Saturated model. m0 <- lmer(y ~ palha * nit * pot + (1 | bloc/ue), data = sugarcane_straw, REML = FALSE) # Simple diagnostic. plot(m0) # qqnorm(residuals(m0, type = "pearson")) # Estimates of the variance components. VarCorr(m0) # Tests for the fixed effects. anova(m0) # Estimates of the effects and fitting measures. # summary(m0) # Final model. mod <- m0 #----------------------------------------------------------------------- # Fitted values. # Experimental values to get estimates. grid <- with(sugarcane_straw, expand.grid(palha = levels(palha), # N = seq(min(N), max(N), length.out = 10), # K = seq(min(N), max(N), length.out = 3), nit = levels(nit), pot = levels(pot), KEEP.OUT.ATTRS = FALSE)) # Matrix of fixed effects. X <- model.matrix(terms(mod), data = cbind(grid, y = 0)) # Confidence intervals. ci <- confint(glht(mod, linfct = X), calpha = univariate_calpha())$confint colnames(ci)[1] <- "fit" grid <- cbind(grid, ci) str(grid) # Sample averages. # averages <- aggregate(nobars(formula(mod)), # data = sugarcane_straw, FUN = mean) # # xyplot(y ~ N | factor(K), groups = palha, data = averages, # type = "o", as.table = TRUE) xyplot(fit ~ pot | nit, groups = palha, data = grid, as.table = TRUE, layout = c(NA, 1), type = "o", xlab = leg$N, ylab = leg$y, auto.key = list(title = leg$palha[1], cex.title = 1.1, text = c(leg$palha[-1]), columns = 2, points = FALSE, lines = TRUE), uy = grid$upr, ly = grid$lwr, desloc = 0.2 * scale(as.integer(grid$palha), scale = FALSE), cty = "bars", length = 0.02, alpha = 0.25, fill = "gray50", prepanel = prepanel.cbH, panel.groups = panel.cbH, panel = panel.superpose)
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"), "----------------------------------------", sep = "\n") sessionInfo()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.