Reproducible Data Analysis of Scientific Cooperations
|
Doses de Nitrogênio e Potássio para Cana-de-açúcar em Diferentes Sistemas de Manejo no Estado do Paraná
2016-06-21
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(wzCoop)Exploratory data analysis
# Object structure.
str(sugarcane_straw)## 'data.frame': 250 obs. of 15 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ bloc : int 1 2 3 4 5 1 2 3 4 5 ...
## $ K : int 0 0 0 0 0 0 0 0 0 0 ...
## $ N : int 0 0 0 0 0 45 45 45 45 45 ...
## $ ncm : num 6.26 6.47 7.98 8.6 7.48 7.14 7.66 7.99 7.31 7.68 ...
## $ pmc : num 1.05 1.22 1.26 1.2 1.2 1.39 1.14 1.26 1.3 1.29 ...
## $ tch : num 54.8 65.5 83.5 86 74.8 ...
## $ pol : num 14.7 15.2 14.9 15.6 16.1 ...
## $ tsh : num 8.02 9.99 12.45 13.45 12.01 ...
## $ tfn : num 17.3 17.3 18.1 18.9 15.8 ...
## $ tfk : num 16.8 18.6 14.7 15.9 15 ...
## $ ue : Factor w/ 10 levels "1.1","2.1","3.1",..: 1 2 3 4 5 1 2 3 4 5 ...
## $ pot : Factor w/ 5 levels "0","45","90",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ nit : Factor w/ 5 levels "0","45","90",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ y : num 16.8 18.6 14.7 15.9 15 ...# Frequencies.
ftable(xtabs(~palha + N + K, data = sugarcane_straw))## K 0 45 90 135 180
## palha N
## 1 0 5 5 5 5 5
## 45 5 5 5 5 5
## 90 5 5 5 5 5
## 135 5 5 5 5 5
## 180 5 5 5 5 5
## 2 0 5 5 5 5 5
## 45 5 5 5 5 5
## 90 5 5 5 5 5
## 135 5 5 5 5 5
## 180 5 5 5 5 5# Checking if is a complete cases dataset.
all(complete.cases(sugarcane_straw))## [1] TRUE# Descriptive measures.
summary(sugarcane_straw)## palha bloc K N ncm
## 1:125 Min. :1 Min. : 0 Min. : 0 Min. :5.150
## 2:125 1st Qu.:2 1st Qu.: 45 1st Qu.: 45 1st Qu.:7.665
## Median :3 Median : 90 Median : 90 Median :7.970
## Mean :3 Mean : 90 Mean : 90 Mean :7.954
## 3rd Qu.:4 3rd Qu.:135 3rd Qu.:135 3rd Qu.:8.287
## Max. :5 Max. :180 Max. :180 Max. :9.720
##
## pmc tch pol tsh
## Min. :0.780 Min. : 50.16 Min. :12.00 Min. : 7.53
## 1st Qu.:1.120 1st Qu.: 72.76 1st Qu.:14.90 1st Qu.:10.89
## Median :1.265 Median : 83.81 Median :15.20 Median :12.93
## Mean :1.273 Mean : 84.62 Mean :15.25 Mean :12.90
## 3rd Qu.:1.438 3rd Qu.: 96.64 3rd Qu.:15.77 3rd Qu.:14.72
## Max. :1.640 Max. :128.36 Max. :17.90 Max. :21.44
##
## tfn tfk ue pot nit
## Min. : 7.65 Min. : 0.77 1.1 : 25 0 :50 0 :50
## 1st Qu.:15.15 1st Qu.:14.37 2.1 : 25 45 :50 45 :50
## Median :16.57 Median :15.62 3.1 : 25 90 :50 90 :50
## Mean :16.41 Mean :15.50 4.1 : 25 135:50 135:50
## 3rd Qu.:18.08 3rd Qu.:17.18 5.1 : 25 180:50 180:50
## Max. :21.33 Max. :20.13 1.2 : 25
## (Other):100
## y
## Min. : 0.77
## 1st Qu.:14.37
## Median :15.62
## Mean :15.50
## 3rd Qu.:17.18
## Max. :20.13
## # 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 0.072606
## bloc (Intercept) 0.085859
## Residual 0.508323# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 0.5592 0.55924 1 5 2.1643 0.2012146
## nit 7.2096 1.80240 4 240 6.9755 2.503e-05 ***
## pot 2.8398 0.70995 4 240 2.7476 0.0290209 *
## palha:nit 5.6340 1.40851 4 240 5.4511 0.0003231 ***
## palha:pot 1.6954 0.42385 4 240 1.6403 0.1647835
## nit:pot 1.6427 0.10267 16 240 0.3973 0.9823795
## palha:nit:pot 1.4598 0.09124 16 240 0.3531 0.9906184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha * nit + pot + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 17 430.48 490.35 -198.24 396.48
## ..1 53 484.60 671.23 -189.30 378.60 17.885 36 0.995anova(m1)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 0.6025 0.60251 1 5 2.1643 0.2012147
## nit 7.2096 1.80240 4 240 6.4745 5.793e-05 ***
## pot 2.8398 0.70995 4 240 2.5503 0.0398982 *
## palha:nit 5.6340 1.40851 4 240 5.0596 0.0006237 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha * N + K + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 8 415.76 443.93 -199.88 399.76
## ..1 53 484.60 671.23 -189.30 378.60 21.162 45 0.9991anova(m2)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 1.6329 1.6329 1 28.678 5.7861 0.02283 *
## N 6.8750 6.8750 1 240.000 24.3611 1.495e-06 ***
## K 2.5461 2.5461 1 240.000 9.0221 0.00295 **
## palha:N 5.3437 5.3437 1 240.000 18.9352 2.002e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2)## Linear mixed model fit by maximum likelihood t-tests use
## Satterthwaite approximations to degrees of freedom [lmerMod]
## Formula: y ~ palha * N + K + (1 | bloc/ue)
## Data: sugarcane_straw
##
## AIC BIC logLik deviance df.resid
## 415.8 443.9 -199.9 399.8 242
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.5240 -0.4284 0.0609 0.4959 2.5747
##
## Random effects:
## Groups Name Variance Std.Dev.
## ue:bloc (Intercept) 0.004319 0.06572
## bloc (Intercept) 0.007372 0.08586
## Residual 0.282211 0.53124
## Number of obs: 250, groups: ue:bloc, 10; bloc, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.428e+00 1.066e-01 5.128e+01 69.663 < 2e-16 ***
## palha2 2.973e-01 1.236e-01 2.868e+01 2.405 0.02283 *
## N 4.903e-03 7.466e-04 2.400e+02 6.567 3.15e-10 ***
## K 1.586e-03 5.279e-04 2.400e+02 3.004 0.00295 **
## palha2:N -4.595e-03 1.056e-03 2.400e+02 -4.351 2.00e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) palha2 N K
## palha2 -0.580
## N -0.630 0.544
## K -0.446 0.000 0.000
## palha2:N 0.446 -0.769 -0.707 0.000# 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)## 'data.frame': 60 obs. of 6 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ N : num 0 0 20 20 40 40 60 60 80 80 ...
## $ K : num 0 0 0 0 0 0 0 0 0 0 ...
## $ fit : num 7.43 7.73 7.53 7.73 7.62 ...
## $ lwr : num 7.22 7.52 7.33 7.54 7.45 ...
## $ upr : num 7.64 7.93 7.72 7.92 7.8 ...# 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 0.0203082
## bloc (Intercept) 0.0052398
## Residual 0.0795324# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 1.96014 1.96014 1 5 309.885 1.085e-05 ***
## nit 1.38310 0.34578 4 240 54.665 < 2.2e-16 ***
## pot 0.58179 0.14545 4 240 22.994 4.441e-16 ***
## palha:nit 0.02894 0.00724 4 240 1.144 0.336548
## palha:pot 0.02131 0.00533 4 240 0.842 0.499621
## nit:pot 0.21837 0.01365 16 240 2.158 0.006996 **
## palha:nit:pot 0.03237 0.00202 16 240 0.320 0.994604
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha + nit * pot + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 29 -475.54 -373.42 266.77 -533.54
## ..1 53 -440.26 -253.62 273.13 -546.26 12.719 24 0.9705anova(m1)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 2.06682 2.06682 1 5 309.885 1.085e-05 ***
## nit 1.38310 0.34577 4 240 51.843 < 2.2e-16 ***
## pot 0.58179 0.14545 4 240 21.807 2.220e-15 ***
## nit:pot 0.21837 0.01365 16 240 2.046 0.01136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha + N * K + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 8 -490.31 -462.14 253.16 -506.31
## ..1 53 -440.26 -253.62 273.13 -546.26 39.948 45 0.6854anova(m2)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 2.31513 2.31513 1 5 309.885 1.085e-05 ***
## N 1.00829 1.00829 1 240 134.961 < 2.2e-16 ***
## K 0.55912 0.55912 1 240 74.839 6.661e-16 ***
## N:K 0.17424 0.17424 1 240 23.322 2.442e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2)## Linear mixed model fit by maximum likelihood t-tests use
## Satterthwaite approximations to degrees of freedom [lmerMod]
## Formula: y ~ palha + N * K + (1 | bloc/ue)
## Data: sugarcane_straw
##
## AIC BIC logLik deviance df.resid
## -490.3 -462.1 253.2 -506.3 242
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.54433 -0.53263 0.04473 0.63706 3.13342
##
## Random effects:
## Groups Name Variance Std.Dev.
## ue:bloc (Intercept) 3.666e-04 0.01915
## bloc (Intercept) 2.746e-05 0.00524
## Residual 7.471e-03 0.08643
## Number of obs: 250, groups: ue:bloc, 10; bloc, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.198e+00 1.943e-02 6.598e+01 61.657 < 2e-16 ***
## palha2 -2.872e-01 1.631e-02 5.000e+00 -17.604 1.08e-05 ***
## N 1.728e-03 1.488e-04 2.400e+02 11.617 < 2e-16 ***
## K 1.287e-03 1.488e-04 2.400e+02 8.651 8.88e-16 ***
## N:K -6.519e-06 1.350e-06 2.400e+02 -4.829 2.44e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) palha2 N K
## palha2 -0.420
## N -0.689 0.000
## K -0.689 0.000 0.667
## N:K 0.563 0.000 -0.816 -0.816
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling# 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)## 'data.frame': 60 obs. of 6 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ N : num 0 0 20 20 40 40 60 60 80 80 ...
## $ K : num 0 0 0 0 0 0 0 0 0 0 ...
## $ fit : num 1.198 0.911 1.233 0.946 1.267 ...
## $ lwr : num 1.16 0.873 1.198 0.911 1.236 ...
## $ upr : num 1.236 0.949 1.267 0.98 1.299 ...# 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 2.4200
## bloc (Intercept) 1.2702
## Residual 7.9438# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 7971.8 7971.8 1 5 126.328 9.737e-05 ***
## nit 11639.3 2909.8 4 240 46.112 < 2.2e-16 ***
## pot 4494.9 1123.7 4 240 17.808 8.122e-13 ***
## palha:nit 494.8 123.7 4 240 1.960 0.1013
## palha:pot 411.0 102.8 4 240 1.628 0.1678
## nit:pot 1195.5 74.7 16 240 1.184 0.2813
## palha:nit:pot 168.3 10.5 16 240 0.167 0.9999
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha + nit + pot + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 13 1818.8 1864.6 -896.40 1792.8
## ..1 53 1865.3 2051.9 -879.65 1759.3 33.514 40 0.7558anova(m1)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 9166.4 9166.4 1 5 126.328 9.737e-05 ***
## nit 11639.3 2909.8 4 240 40.102 < 2.2e-16 ***
## pot 4494.9 1123.7 4 240 15.487 2.768e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha + N + K + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 7 1820.2 1844.9 -903.11 1806.2
## ..1 53 1865.3 2051.9 -879.65 1759.3 46.932 46 0.4341anova(m2)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 9693.5 9693.5 1 5 126.33 9.737e-05 ***
## N 11187.3 11187.3 1 240 145.80 < 2.2e-16 ***
## K 3945.6 3945.6 1 240 51.42 9.186e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2)## Linear mixed model fit by maximum likelihood t-tests use
## Satterthwaite approximations to degrees of freedom [lmerMod]
## Formula: y ~ palha + N + K + (1 | bloc/ue)
## Data: sugarcane_straw
##
## AIC BIC logLik deviance df.resid
## 1820.2 1844.9 -903.1 1806.2 243
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.9369 -0.5248 0.0271 0.6189 2.6764
##
## Random effects:
## Groups Name Variance Std.Dev.
## ue:bloc (Intercept) 5.311 2.305
## bloc (Intercept) 1.613 1.270
## Residual 76.733 8.760
## Number of obs: 250, groups: ue:bloc, 10; bloc, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 79.830440 1.796238 25.010000 44.443 < 2e-16 ***
## palha2 -20.578400 1.830887 5.000000 -11.240 9.74e-05 ***
## N 0.105115 0.008705 240.000000 12.075 < 2e-16 ***
## K 0.062425 0.008705 240.000000 7.171 9.19e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) palha2 N
## palha2 -0.510
## N -0.436 0.000
## K -0.436 0.000 0.000# 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)## 'data.frame': 60 obs. of 6 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ N : num 0 0 20 20 40 40 60 60 80 80 ...
## $ K : num 0 0 0 0 0 0 0 0 0 0 ...
## $ fit : num 79.8 59.3 81.9 61.4 84 ...
## $ lwr : num 76.3 55.7 78.5 58 80.8 ...
## $ upr : num 83.4 62.8 85.3 64.7 87.3 ...# 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 0.59515
## bloc (Intercept) 0.00000
## Residual 0.92467# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 0.000261 0.0002605 1 10 0.0003047 0.9864
## nit 0.053042 0.0132606 4 240 0.0155092 0.9995
## pot 0.038214 0.0095536 4 240 0.0111736 0.9998
## palha:nit 0.003798 0.0009494 4 240 0.0011104 1.0000
## palha:pot 0.036210 0.0090524 4 240 0.0105874 0.9998
## nit:pot 0.070650 0.0044156 16 240 0.0051644 1.0000
## palha:nit:pot 0.079654 0.0049784 16 240 0.0058226 1.0000# 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 0.68218
## bloc (Intercept) 0.00000
## Residual 1.37393# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 86.426 86.426 1 10 45.784 4.941e-05 ***
## nit 278.180 69.545 4 240 36.841 < 2.2e-16 ***
## pot 103.586 25.896 4 240 13.719 4.334e-10 ***
## palha:nit 10.766 2.692 4 240 1.426 0.2260
## palha:pot 8.375 2.094 4 240 1.109 0.3528
## nit:pot 29.942 1.871 16 240 0.991 0.4668
## palha:nit:pot 3.930 0.246 16 240 0.130 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha + nit + pot + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 13 940.56 986.33 -457.28 914.56
## ..1 53 994.00 1180.63 -444.00 888.00 26.558 40 0.9492anova(m1)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 96.539 96.539 1 10 45.784 4.941e-05 ***
## nit 278.180 69.545 4 240 32.982 < 2.2e-16 ***
## pot 103.586 25.896 4 240 12.281 4.218e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## Data: sugarcane_straw
## Models:
## object: y ~ palha + N + K + (1 | bloc/ue)
## ..1: y ~ palha * nit * pot + (1 | bloc/ue)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## object 7 939.43 964.08 -462.71 925.43
## ..1 53 994.00 1180.63 -444.00 888.00 37.429 46 0.812anova(m2)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 101.012 101.012 1 10 45.784 4.941e-05 ***
## N 268.337 268.337 1 240 121.624 < 2.2e-16 ***
## K 89.981 89.981 1 240 40.784 8.729e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m2)## Linear mixed model fit by maximum likelihood t-tests use
## Satterthwaite approximations to degrees of freedom [lmerMod]
## Formula: y ~ palha + N + K + (1 | bloc/ue)
## Data: sugarcane_straw
##
## AIC BIC logLik deviance df.resid
## 939.4 964.1 -462.7 925.4 243
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2878 -0.6152 0.0526 0.5537 2.6531
##
## Random effects:
## Groups Name Variance Std.Dev.
## ue:bloc (Intercept) 0.4526 0.6728
## bloc (Intercept) 0.0000 0.0000
## Residual 2.2063 1.4854
## Number of obs: 250, groups: ue:bloc, 10; bloc, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 12.160400 0.378781 17.510000 32.104 < 2e-16 ***
## palha2 -3.147280 0.465134 10.000000 -6.766 4.94e-05 ***
## N 0.016280 0.001476 240.000000 11.028 < 2e-16 ***
## K 0.009427 0.001476 240.000000 6.386 8.73e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) palha2 N
## palha2 -0.614
## N -0.351 0.000
## K -0.351 0.000 0.000# 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)## 'data.frame': 60 obs. of 6 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ N : num 0 0 20 20 40 40 60 60 80 80 ...
## $ K : num 0 0 0 0 0 0 0 0 0 0 ...
## $ fit : num 12.16 9.01 12.49 9.34 12.81 ...
## $ lwr : num 11.42 8.27 11.76 8.61 12.1 ...
## $ upr : num 12.9 9.76 13.21 10.06 13.52 ...# 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 0.31359
## bloc (Intercept) 0.17175
## Residual 1.43422# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 46.72 46.718 1 5 22.7118 0.0050340 **
## nit 48.12 12.031 4 240 5.8489 0.0001656 ***
## pot 26.94 6.734 4 240 3.2739 0.0122767 *
## palha:nit 75.22 18.804 4 240 9.1417 6.860e-07 ***
## palha:pot 24.78 6.194 4 240 3.0114 0.0188883 *
## nit:pot 370.94 23.184 16 240 11.2707 < 2.2e-16 ***
## palha:nit:pot 215.09 13.443 16 240 6.5355 3.458e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## 'data.frame': 50 obs. of 6 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ nit : Factor w/ 5 levels "0","45","90",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ pot : Factor w/ 5 levels "0","45","90",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ fit : num 17.5 17.2 19.3 15.5 13.4 ...
## $ lwr : num 16.2 15.9 18 14.3 12.1 ...
## $ upr : num 18.8 18.5 20.5 16.8 14.7 ...# 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)## Groups Name Std.Dev.
## ue:bloc (Intercept) 0.26525
## bloc (Intercept) 0.15615
## Residual 1.15447# Tests for the fixed effects.
anova(m0)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## palha 160.385 160.385 1 5 120.338 0.0001095 ***
## nit 23.789 5.947 4 240 4.462 0.0016995 **
## pot 139.954 34.988 4 240 26.252 < 2.2e-16 ***
## palha:nit 9.992 2.498 4 240 1.874 0.1155932
## palha:pot 71.767 17.942 4 240 13.462 6.49e-10 ***
## nit:pot 281.636 17.602 16 240 13.207 < 2.2e-16 ***
## palha:nit:pot 236.242 14.765 16 240 11.078 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 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)## 'data.frame': 50 obs. of 6 variables:
## $ palha: Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ nit : Factor w/ 5 levels "0","45","90",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ pot : Factor w/ 5 levels "0","45","90",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ fit : num 16.2 14 15.9 14.5 18 ...
## $ lwr : num 15.1 12.9 14.9 13.4 17 ...
## $ upr : num 17.2 15 17 15.5 19 ...# 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
## terça, 21 de junho de 2016, 22:40
## ----------------------------------------
## R version 3.3.0 (2016-05-03)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
##
## locale:
## [1] LC_CTYPE=en_US.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] stats graphics grDevices utils datasets methods
## [7] base
##
## other attached packages:
## [1] wzCoop_0.0-1 wzRfun_0.65 latticeExtra_0.6-28
## [4] RColorBrewer_1.1-2 knitr_1.12.3 multcomp_1.4-5
## [7] TH.data_1.0-7 MASS_7.3-45 survival_2.39-4
## [10] mvtnorm_1.0-5 doBy_4.5-15 lmerTest_2.0-30
## [13] lme4_1.1-12 Matrix_1.2-6 gridExtra_2.2.1
## [16] lattice_0.20-33 devtools_1.11.1
##
## loaded via a namespace (and not attached):
## [1] zoo_1.7-13 splines_3.3.0 colorspace_1.2-6
## [4] testthat_1.0.2 htmltools_0.3.5 rpanel_1.1-3
## [7] yaml_2.1.13 chron_2.3-47 nloptr_1.0.4
## [10] foreign_0.8-66 withr_1.0.1 plyr_1.8.4
## [13] stringr_1.0.0 munsell_0.4.3 gtable_0.2.0
## [16] codetools_0.2-14 memoise_1.0.0 evaluate_0.9
## [19] curl_0.9.7 Rcpp_0.12.4 acepack_1.3-3.3
## [22] scales_0.4.0 formatR_1.3 Hmisc_3.17-4
## [25] ggplot2_2.1.0 digest_0.6.9 stringi_1.0-1
## [28] grid_3.3.0 tools_3.3.0 sandwich_2.3-4
## [31] magrittr_1.5 Formula_1.2-1 cluster_2.0.4
## [34] crayon_1.3.1 data.table_1.9.6 minqa_1.2.4
## [37] rmarkdown_0.9.6 roxygen2_5.0.1 httr_1.1.0
## [40] rpart_4.1-10 R6_2.1.2 nnet_7.3-12
## [43] nlme_3.1-128 git2r_0.14.0
wzCoop - Reproducible Data Analysis of Scientific
Cooperations
|
Walmes Zeviani |