Definições da sessão
#-----------------------------------------------------------------------
# Prof. Dr. Walmes M. Zeviani
# leg.ufpr.br/~walmes · github.com/walmes
# walmes@ufpr.br · @walmeszeviani
# Laboratory of Statistics and Geoinformation (LEG)
# Department of Statistics · Federal University of Paraná
# 2020-nov-16 · Curitiba/PR/Brazil
#-----------------------------------------------------------------------
#-----------------------------------------------------------------------
# Pacotes.
rm(list = objects())
library(nlme)
library(multcomp)
library(emmeans)
library(tidyverse)
library(readxl)
# source("https://raw.githubusercontent.com/walmes/wzRfun/master/R/pairwise.R")
Importação e análise exploratória
#-----------------------------------------------------------------------
# Importação e preparo.
tb <- read_xlsx("Volatilização_v1.xlsx", sheet = "volatilizacao")
str(tb)
## tibble [90 × 9] (S3: tbl_df/tbl/data.frame)
## $ tratamento: num [1:90] 1 1 1 2 2 2 3 3 3 4 ...
## $ bloco : num [1:90] 1 2 3 1 2 3 1 2 3 1 ...
## $ tipo_ureia: chr [1:90] "ureia" "ureia" "ureia" "ureia" ...
## $ tipo_mo : chr [1:90] "ma" "ma" "ma" "ma" ...
## $ dose_mo : num [1:90] 0 0 0 150 150 150 300 300 300 600 ...
## $ perd_ac_% : num [1:90] 53.3 51.5 48.2 51.8 48.4 ...
## $ ph : num [1:90] 8.05 7.84 8.5 7.94 7.71 7.74 8.02 7.78 7.6 7.74 ...
## $ ph_ca : num [1:90] 7.44 7.33 7.18 7.21 7.06 7.22 7.31 7.03 7.07 7.33 ...
## $ ph_agua : num [1:90] 7.33 7.25 7.14 7.12 7.03 7.14 7.24 7 7.05 7.21 ...
tb <- tb %>%
mutate_at(c("tratamento", "bloco", "tipo_ureia", "tipo_mo"),
factor) %>%
rename("nitro" = 6) %>%
rename_at(c("tipo_ureia", "tipo_mo"), str_to_title) %>%
mutate(Dose_mo = factor(dose_mo))
tb_long <- tb %>%
gather(key = "metodo", value = "ph", starts_with("ph"))
str(tb_long)
## tibble [270 × 9] (S3: tbl_df/tbl/data.frame)
## $ tratamento: Factor w/ 30 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ bloco : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
## $ Tipo_ureia: Factor w/ 3 levels "incorporado",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Tipo_mo : Factor w/ 2 levels "ma","tmo": 1 1 1 1 1 1 1 1 1 1 ...
## $ dose_mo : num [1:270] 0 0 0 150 150 150 300 300 300 600 ...
## $ nitro : num [1:270] 53.3 51.5 48.2 51.8 48.4 ...
## $ Dose_mo : Factor w/ 5 levels "0","150","300",..: 1 1 1 2 2 2 3 3 3 4 ...
## $ metodo : chr [1:270] "ph" "ph" "ph" "ph" ...
## $ ph : num [1:270] 8.05 7.84 8.5 7.94 7.71 7.74 8.02 7.78 7.6 7.74 ...
#-----------------------------------------------------------------------
# Análise exploratória.
# Relação nos dados crus.
ggplot(data = tb_long,
mapping = aes(x = ph, y = nitro)) +
facet_wrap(facets = ~metodo, ncol = 1) +
# geom_point(mapping = aes(color = tipo_mo)) +
# geom_point(mapping = aes(color = tipo_ureia)) +
# geom_point(mapping = aes(color = dose_mo)) +
geom_point() +
geom_smooth(method = "lm") +
ggtitle("Análise de correlação com as respostas")
## `geom_smooth()` using formula 'y ~ x'
Análise de correlação
# ATTENTION: qual a relação após descontar o efeito dos termos
# experimentais?
# Declara modelo de MANOVA com preditor máximo pra estudar a correlação
# nos resíduos livre dos eventuais efeitos de termos impostos pelo
# experimento.
m0 <- lm(cbind(nitro, ph, ph_ca, ph_agua) ~
bloco + Tipo_ureia * Tipo_mo * Dose_mo,
data = tb)
# Dispersão dos resíduos.
lattice::splom(residuals(m0), type = c("p", "r"))
# Correlação de Pearson.
rc <- Hmisc::rcorr(residuals(m0))
rc
## nitro ph ph_ca ph_agua
## nitro 1.00 0.38 -0.17 -0.18
## ph 0.38 1.00 -0.26 -0.23
## ph_ca -0.17 -0.26 1.00 0.78
## ph_agua -0.18 -0.23 0.78 1.00
##
## n= 90
##
##
## P
## nitro ph ph_ca ph_agua
## nitro 0.0003 0.1038 0.0877
## ph 0.0003 0.0143 0.0271
## ph_ca 0.1038 0.0143 0.0000
## ph_agua 0.0877 0.0271 0.0000
M <- rc$r
p_mat <- rc$P
# Correlograma.
corrplot::corrplot(M,
method = "ellipse",
addCoef.col = "black",
number.cex = 0.75,
addCoefasPercent = FALSE,
p.mat = p_mat,
sig.level = 0.01)
rect(xleft = rep(0.5, 2),
ybottom = c(0.5, ncol(M) - 0.5),
xright = c(1.5, ncol(M) + 0.5),
ytop = rep(ncol(M) + 0.5, 2),
col = NA,
border = "red",
lwd = 2)
res_long <- residuals(m0) %>%
as.data.frame() %>%
gather(key = "metodo", valu = "ph", -1)
# Análise de correlação de Pearson com os resíduos.
tb_cor <- res_long %>%
group_by(metodo) %>%
do(broom::tidy(cor.test(.$nitro, .$ph))) %>%
ungroup() %>%
select(-method, -alternative, -parameter) %>%
mutate(sig = cut(p.value,
c(0, 0.001, 0.01, 0.05, 0.1, 1),
c("***", "**", "*", "·", "")))
tb_cor
## # A tibble: 3 x 7
## metodo estimate statistic p.value conf.low conf.high sig
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
## 1 ph 0.376 3.81 0.000257 0.183 0.541 "***"
## 2 ph_agua -0.181 -1.73 0.0877 -0.374 0.0271 "·"
## 3 ph_ca -0.173 -1.64 0.104 -0.367 0.0358 ""
ggplot(data = res_long,
mapping = aes(x = ph, y = nitro)) +
facet_wrap(facets = ~metodo, ncol = 1) +
geom_vline(xintercept = 0, lty = 2, size = 0.25) +
geom_hline(yintercept = 0, lty = 2, size = 0.25) +
geom_point() +
geom_label(data = tb_cor,
mapping = aes(x = 0.25,
y = -6,
label = sprintf("r = %0.2f (%0.4f%s)",
estimate,
p.value, sig))) +
geom_smooth(method = "lm", col = "orange") +
ggtitle("Análise de correlação com os resíduos") +
labs(x = "pH do solo determinado por algum método",
y = "Nitrogênio liberado acumulado",
caption = "Valores de N e pH observados na última coleta.")
## `geom_smooth()` using formula 'y ~ x'
