In walmes/RDASC: Reproducible Data Analysis of Cientific Cooperations
Session Definition
# 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")
Exploratory data analysis
# 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)
})
Número de colmos por metro
#-----------------------------------------------------------------------
# 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)
Peso médio de colmo
#-----------------------------------------------------------------------
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)
Produção de cana-de-açúcar
#-----------------------------------------------------------------------
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)
Teor de sacarose aparente
#-----------------------------------------------------------------------
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)
Produção de sacarose
#-----------------------------------------------------------------------
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)
Teor de nitrogênio nas folhas
#-----------------------------------------------------------------------
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.
Teor de potássio nas folhas
#-----------------------------------------------------------------------
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)
Session information
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"),
"----------------------------------------", sep = "\n")
sessionInfo()
walmes/RDASC documentation built on Jan. 10, 2021, 8:02 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Session Definition
# 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")
Exploratory data analysis
# 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) })
Número de colmos por metro
#----------------------------------------------------------------------- # 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)
Peso médio de colmo
#----------------------------------------------------------------------- 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)
Produção de cana-de-açúcar
#----------------------------------------------------------------------- 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)
Teor de sacarose aparente
#----------------------------------------------------------------------- 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)
Produção de sacarose
#----------------------------------------------------------------------- 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)
Teor de nitrogênio nas folhas
#----------------------------------------------------------------------- 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.
Teor de potássio nas folhas
#----------------------------------------------------------------------- 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)
Session information
cat(format(Sys.time(), format = "%A, %d de %B de %Y, %H:%M"), "----------------------------------------", sep = "\n") sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.