Ajuste de curvas de liberação de nitrogênio com duas fases
Lílian Angélica Moreira & Walmes Marques Zeviani
2021-03-09
1 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)2 Importação e análise exploratória
#-----------------------------------------------------------------------
# Importação e preparo.
# Leitura.
tb <- read_xlsx("Volatilização_v1.xlsx", sheet = "N perd. acumulado (%)")
# str(tb)
# Conversão de variáveis para fator.
tb <- tb %>%
mutate_at(c("tratamento", "bloco", "tipo_ureia", "tipo_mo"), factor)
# Empilha as colunas dos dias e define a unidade experimental.
tb <- tb %>%
gather(key = "dia", value = "nitro", -(tratamento:dose_mo)) %>%
mutate(dia = as.integer(dia),
ue = interaction(tratamento, bloco))
str(tb)## tibble [4,140 × 8] (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:4140] 0 0 0 150 150 150 300 300 300 600 ...
## $ dia : int [1:4140] 1 1 1 1 1 1 1 1 1 1 ...
## $ nitro : num [1:4140] 0.0809 0.0644 0.0364 0.1191 0.0448 ...
## $ ue : Factor w/ 90 levels "1.1","2.1","3.1",..: 1 31 61 2 32 62 3 33 63 4 ...#-----------------------------------------------------------------------
# Análise exploratória.
ggplot(data = tb,
mapping = aes(x = dia, y = nitro, color = tipo_mo)) +
facet_grid(facets = dose_mo ~ tipo_ureia) +
geom_point() +
# scale_x_sqrt() +
# scale_x_log10() +
stat_summary(mapping = aes(group = tipo_mo),
geom = "line", fun = "mean")
3 Modelo para a fase de pronta liberação
- https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)
- https://www.mdpi.com/2073-8994/11/2/204/htm
- https://www.researchgate.net/publication/257417959_The_Hill_equation_and_the_origin_of_quantitative_pharmacology
#-----------------------------------------------------------------------
# Modelo Hill-Langmuir (HL).
#
# * Intercepto é 0.
# * É um modelo sigmoidal.
# * Parâmetro de assíntota e meia-vida.
# * É possível determinar a taxa máxima.
parlist <- list(tha = 25, thv = 6, thk = 2.5)
with(parlist, {
# Equação 1.
curve(tha * x^thk/(thv^thk + x^thk),
from = 0, to = 30, col = "cyan", lwd = 3)
# Equação 2.
curve(tha/(1 + (thv/x)^thk), add = TRUE,
from = 0, to = 30, col = "orange",
lwd = 3, lty = 2)
abline(v = thv, h = tha/2, lty = 2, col = "gray")
})
# TODO:
# Qual a taxa máxima do HL?
# Derivada de primeira ordem.
D(expression(tha/(1 + (thv/x)^thk)), "x")## tha * ((thv/x)^(thk - 1) * (thk * (thv/x^2)))/(1 + (thv/x)^thk)^2with(parlist, {
curve(tha * ((thv/x)^(thk - 1) * (thk * (thv/x^2)))/(1 + (thv/x)^thk)^2,
from = 0, to = 30, col = "cyan", lwd = 3, ylab = "f'(x)")
# abline(v = thv, h = tha/2, lty = 2, col = "gray")
})
# Derivada de segunda ordem.
D(D(expression(tha/(1 + (thv/x)^thk)), "x"), "x")## -(tha * ((thv/x)^(thk - 1) * (thk * (thv * (2 * x)/(x^2)^2)) +
## (thv/x)^((thk - 1) - 1) * ((thk - 1) * (thv/x^2)) * (thk *
## (thv/x^2)))/(1 + (thv/x)^thk)^2 - tha * ((thv/x)^(thk -
## 1) * (thk * (thv/x^2))) * (2 * ((thv/x)^(thk - 1) * (thk *
## (thv/x^2)) * (1 + (thv/x)^thk)))/((1 + (thv/x)^thk)^2)^2)with(parlist, {
curve(-(tha * ((thv/x)^(thk - 1) * (thk * (thv * (2 * x)/(x^2)^2)) +
(thv/x)^((thk - 1) - 1) * ((thk - 1) * (thv/x^2)) * (thk *
(thv/x^2)))/(1 + (thv/x)^thk)^2 - tha * ((thv/x)^(thk -
1) * (thk * (thv/x^2))) * (2 * ((thv/x)^(thk - 1) * (thk *
(thv/x^2)) * (1 + (thv/x)^thk)))/((1 + (thv/x)^thk)^2)^2),
from = 0, to = 30, col = "cyan", lwd = 3, ylab = "f''(x)")
abline(v = thv, h = tha/2, lty = 2, col = "gray")
abline(h = 0, col = "red")
})
# Usando o Wolfram Alpha com o seguinte argumento de consulta
# solve D[1/(1 + (v/x)^k), {x, 2}] = 0 wrt x
#
# O resultado é
# x = e^(-(i π)/k) (-k^2 - k)^(-1/k) (k^2 - k)^(1/k) v
# ATTENTION: Para determinar o ponto de inflexão e a taxa máxima, vamos
# usar métodos númericos pois os resultados serão os mesmos.Dessa forma o modelo a ser usado para a \[ f(t) = \frac{\theta_a}{1 + \left(\frac{\theta_v}{t}\right)^{\theta_k}}. \]
em que
- \(\theta_a\) é o conteúdo total de fácil liberação.
- \(\theta_v\) é a EC50, ou seja, ponto no domínio de \(t\) que corresponde à 50% de \(\theta_a\).
- \(\theta_k\) é um parâmetro de forma associado à taxa máxima de liberação.
4 Modelo linear para a fase de lenta liberação
O modelo ajustado é uma soma de duas funções. A função Hill-Langmuir é usada para descrever o comportamento do nitrogênio de fácil ou pronta liberação enquanto que o termo adicional é usado para caracterízar a fase de lenta liberação. A fase de lenta liberação é uma função de taxa constante obtida pela inclusão de um termo linear. Para que haja contínuidade, usa-se um termo quadrático que se solda ao termo linear.
\[ \begin{aligned} f(t) &= f_1 (t) + f_2 (t) \\ f_1(t) &= \frac{\theta_a}{ 1 + \left(\frac{\theta_v}{t} \right)^{\theta_k}} \\ f_2(t) &= \begin{cases} 0 & t, < \theta_b \\ \theta_x (t - \theta_b)^2/2, & \theta_b < t < \theta_b + 1 \\ \theta_x (t - \theta_b) & t > \theta_b + 1 \end{cases}\\ \end{aligned} \]
em que
- \(\theta_a\) é o conteúdo total de fácil liberação.
- \(\theta_v\) é a EC_50, ou seja, ponto no domínio de \(t\) que corresponde à 50% de \(\theta_a\).
- \(\theta_k\) é um parâmetro de forma associado à taxa máxima de liberação.
- \(\theta_b\) é o ponto em que inicia a liberação da fase de difícil liberação (break).
- \(\theta_x\) é a taxa da fase de difícil liberação.
#-----------------------------------------------------------------------
# Gráfico da função para a fase de lenta liberação.
with(list(thx = 2, thb = 1), {
curve(0 +
(thx/2 * (x - thb)^2) * (x >= thb) * (x < (thb + 1))+
(thx/2 + thx * (x - (thb + 1))) * (x >= (thb + 1)),
from = 0, to = 4, col = "gray")
curve(0 * x, to = thb, col = "orange",
lwd = 2, lty = 2, add = TRUE)
curve((thx/2 * (x - thb)^2),
from = thb, to = thb + 1,
col = "magenta",
lwd = 2, lty = 2, add = TRUE)
curve((thx/2+ thx * (x - (thb + 1))),
from = thb + 1,
col = "cyan",
lwd = 2, lty = 2, add = TRUE)
abline(v = thb + 0:1, lty = 2, col = "red")
abline(h = thx/2, lty = 2, col = "red")
})
5 Ajuste para uma unidade experimental (opcional)
#-----------------------------------------------------------------------
# Aplicação em uma unidade experimental.
# Testa os valores iniciais.
i <- "4"
tbi <- tb %>%
filter(tratamento == i, bloco == "1", dia <= 30) %>%
arrange(dia)
start <- list(tha = 30, thv = 7, thk = 2.4, thx = 0.8, thb = 12)
model <- nitro ~ tha/(1 + (thv/dia)^thk) +
(thx/2 * (dia - thb)^2) * (dia >= thb) * (dia < (thb + 1))+
(thx/2 + thx * (dia - (thb + 1))) * (dia >= (thb + 1))
plot(nitro ~ dia, data = tbi)
with(start, {
curve(tha/(1 + (thv/dia)^thk) +
(thx/2 * (dia - thb)^2) * (dia >= thb) * (dia < (thb + 1))+
(thx/2 + thx * (dia - (thb + 1))) * (dia >= (thb + 1)),
xname = "dia",
add = TRUE,
col = "red")
curve(tha/(1 + (thv/x)^thk),
col = "orange",
lwd = 2, lty = 2, add = TRUE)
curve((thx/2 * (x - thb)^2),
from = thb, to = thb + 1,
col = "magenta",
lwd = 2, lty = 2, add = TRUE)
curve((thx/2+ thx * (x - (thb + 1))),
from = thb + 1,
col = "cyan",
lwd = 2, lty = 2, add = TRUE)
})
# Ajuste do modelo.
n0 <- nls(model,
data = tbi,
start = start)
summary(n0)##
## Formula: nitro ~ tha/(1 + (thv/dia)^thk) + (thx/2 * (dia - thb)^2) * (dia >=
## thb) * (dia < (thb + 1)) + (thx/2 + thx * (dia - (thb + 1))) *
## (dia >= (thb + 1))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## tha 36.40758 0.74766 48.70 <2e-16 ***
## thv 7.76876 0.15185 51.16 <2e-16 ***
## thk 2.70078 0.08654 31.21 <2e-16 ***
## thx 0.67927 0.03239 20.97 <2e-16 ***
## thb 14.94997 0.50415 29.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3439 on 25 degrees of freedom
##
## Number of iterations to convergence: 8
## Achieved convergence tolerance: 8.467e-06# Confere o ajuste.
plot(nitro ~ dia, data = tbi)
with(as.list(coef(n0)), {
curve(tha/(1 + (thv/dia)^thk) +
(thx/2 * (dia - thb)^2) * (dia >= thb) * (dia < (thb + 1))+
(thx/2 + thx * (dia - (thb + 1))) * (dia >= (thb + 1)),
xname = "dia",
add = TRUE,
col = "red")
curve(tha/(1 + (thv/x)^thk),
col = "orange",
lwd = 2, lty = 2, add = TRUE)
curve((thx/2 * (x - thb)^2),
from = thb, to = thb + 1,
col = "magenta",
lwd = 2, lty = 2, add = TRUE)
curve((thx/2+ thx * (x - (thb + 1))),
from = thb + 1,
col = "cyan",
lwd = 2, lty = 2, add = TRUE)
})
# Verifica necessidade de estrutura de correlação.
par(mfrow = c(1, 2))
acf(residuals(n0))
pacf(residuals(n0))
layout(1)
# Ajusta modelo com estrutura de correlação.
n1 <- gnls(model = model,
data = tbi,
start = coef(n0),
correlation = corCAR1(value = 0.5))
# Compara os ajustes.
anova(n1, n0)## Model df AIC BIC logLik Test L.Ratio p-value
## n1 1 7 26.46092 36.26930 -6.230460
## n0 2 6 27.61473 36.02192 -7.807366 1 vs 2 3.153813 0.07576 Ajuste do modelo para cada unidade experimental
#-----------------------------------------------------------------------
# Pela análise exploratória, o modelo terá de ser composto por duas
# funções: uma sigmoidal e uma linear ou suplinear.
# Hill-Langmuir.
start <- list(tha = 34, thv = 9, thk = 3, thx = 0.7, thb = 12)
model <- nitro ~ tha/(1 + (thv/dia)^thk) +
(thx/2 * (dia - thb)^2) * (dia >= thb) * (dia < (thb + 1))+
(thx/2 + thx * (dia - (thb + 1))) * (dia >= (thb + 1)) | ue
# Usar uma janela menor no tempo.
tb_ref <- tb %>%
filter(dia <= 25)
# Ajuste em batelada.
n0 <- nlsList(model = model,
start = start,
data = tb_ref,
control = list(maxiter = 500))
# dput(round(colMeans(coef(n0), na.rm = TRUE), 3))
# summary(n0)
coef(n0)## tha thv thk thx thb
## 1.1 33.58612 6.656874 3.277294 0.7772601 13.065246
## 2.1 35.67208 8.524628 3.330366 0.6817638 15.818588
## 3.1 35.32390 9.236118 3.049976 0.5926979 16.619164
## 4.1 36.52014 7.791831 2.687705 0.6502496 14.835952
## 5.1 33.29081 8.537017 3.374067 0.8514720 15.680890
## 6.1 33.62386 11.385044 3.734068 0.4415414 17.951950
## 7.1 39.77972 10.038958 2.819150 0.5411785 17.581950
## 8.1 33.60321 10.079917 2.971285 0.5644352 17.383849
## 9.1 33.36010 9.552995 3.045847 0.7589567 15.885254
## 10.1 31.89066 9.278894 3.621736 0.6771097 16.576982
## 11.1 22.90775 14.067566 2.903010 0.5670568 11.215957
## 12.1 42.58537 9.798676 2.532839 1.0339008 17.614671
## 13.1 34.98163 7.529597 2.794650 0.7683214 13.254498
## 14.1 33.08408 9.005693 3.015774 0.7450452 15.024359
## 15.1 37.99510 9.690620 3.178435 0.6542898 16.907501
## 16.1 32.52767 8.494935 3.088053 0.7270530 15.892628
## 17.1 34.83127 8.430019 2.581026 0.6910378 15.083291
## 18.1 36.63228 9.861590 3.251149 0.9148590 18.060759
## 19.1 31.28040 7.897410 3.082177 0.7795188 14.272777
## 20.1 32.21535 7.084094 4.147288 0.7940039 12.964548
## 21.1 21.70461 8.993761 3.236277 1.2162000 8.012687
## 22.1 32.67593 7.525220 2.857635 0.8114159 13.433723
## 23.1 28.22322 7.648118 3.819489 1.2085657 13.548602
## 24.1 19.41354 7.363176 3.171855 1.1531235 6.031353
## 25.1 34.16448 6.791822 3.983707 1.0665280 10.280729
## 26.1 18.22907 12.571137 3.904168 0.5994767 9.985410
## 27.1 17.15099 7.851435 3.197153 0.8596073 12.887596
## 28.1 36.38514 9.622932 3.045936 0.5194406 9.370984
## 29.1 28.71565 8.406597 2.884131 0.7276628 6.775502
## 30.1 35.50892 9.297235 3.103074 0.6437141 16.876507
## 1.2 34.73246 9.629716 3.208510 0.7076047 17.001486
## 2.2 33.91661 10.606596 3.956057 0.6363001 17.895668
## 3.2 32.12648 9.931056 3.462516 0.9020426 17.559058
## 4.2 29.37238 8.258710 3.884784 0.8772644 14.431280
## 5.2 39.38964 8.117107 2.984024 0.9402306 12.378662
## 6.2 20.45910 8.697471 3.682591 1.0140742 8.016702
## 7.2 33.19944 10.250488 3.255761 0.6406763 17.301233
## 8.2 24.42983 11.655670 2.411659 0.7194367 13.093323
## 9.2 32.19837 7.818715 3.247134 0.8004738 14.085217
## 10.2 33.39755 9.251151 3.337947 0.7121602 16.104641
## 11.2 20.05289 13.820364 3.257913 0.6851702 11.103964
## 12.2 31.97060 9.923016 3.844516 0.7828333 16.081560
## 13.2 39.13703 9.677983 3.303349 1.0969904 17.584017
## 14.2 34.77704 7.879172 3.038525 0.7863784 14.015503
## 15.2 37.08099 9.737861 3.365233 0.7686973 16.036024
## 16.2 32.30347 10.152834 4.016175 1.0157205 17.269321
## 17.2 33.65243 8.240565 3.465759 0.7864580 14.831451
## 18.2 32.41613 7.690344 2.967316 0.7872858 13.353307
## 19.2 33.20645 8.340270 3.284952 0.7934238 14.848524
## 20.2 32.61874 9.593924 3.241645 0.6123764 16.860560
## 21.2 33.10073 12.975656 2.858706 0.4354943 11.411123
## 22.2 34.69191 7.498969 3.063129 0.8368207 13.867480
## 23.2 32.37690 8.450881 3.540755 0.7867101 14.589495
## 24.2 40.91890 8.961858 2.739989 0.5231081 16.823998
## 25.2 35.76743 9.267451 2.996293 0.6572024 16.688869
## 26.2 21.01858 11.492817 4.194636 0.9622516 9.929209
## 27.2 34.53084 8.888344 3.613617 0.7382474 15.638967
## 28.2 35.33349 8.577633 3.224586 0.7205495 15.549606
## 29.2 27.79830 9.794805 3.419327 0.7062403 8.638370
## 30.2 31.46817 7.435367 3.648534 0.8175070 13.088406
## 1.3 30.83072 10.169361 3.430714 0.9987553 18.229070
## 2.3 39.49564 9.109129 3.299951 0.6936826 17.146501
## 3.3 30.29853 10.960110 3.232348 0.3514892 7.442767
## 4.3 37.24515 8.456981 3.309732 0.7279053 15.889339
## 5.3 34.39000 9.554289 3.005215 0.6171424 16.927811
## 6.3 38.13192 8.310038 3.700873 0.7159167 15.318560
## 7.3 34.56246 9.203107 3.556281 0.6848879 16.577966
## 8.3 32.84542 9.059816 2.985039 0.6951514 15.900856
## 9.3 35.37270 9.216936 3.109625 0.6598567 16.775335
## 10.3 35.52058 8.654899 2.894876 0.7089587 15.872334
## 11.3 15.00973 9.984967 3.835612 0.8387281 10.379424
## 12.3 29.83275 6.066689 3.215025 0.9040997 10.461604
## 13.3 33.38869 8.473070 3.197872 0.7955296 14.853813
## 14.3 34.05190 7.691122 3.257190 0.8894890 13.168346
## 15.3 30.78259 7.406020 3.007076 0.8919132 13.901143
## 16.3 35.55548 7.740155 3.345940 0.8624778 14.127264
## 17.3 37.52716 9.886593 2.649968 0.5424480 17.052858
## 18.3 37.34082 9.949100 2.777021 0.6174814 17.510616
## 19.3 32.43189 8.446989 3.134141 0.7408944 14.993836
## 20.3 29.66239 7.842758 3.703055 0.8030714 13.855119
## 21.3 22.27025 10.215519 3.407450 0.7037748 8.809124
## 22.3 30.52760 8.947921 3.492180 0.7073947 15.154828
## 23.3 34.61722 6.913988 3.746776 0.8798018 12.714240
## 24.3 34.08141 8.819427 3.052510 0.6935386 16.124478
## 25.3 34.18364 8.295390 3.021885 0.7286689 15.034368
## 26.3 30.32916 13.061027 3.202869 0.4289752 18.282201
## 27.3 32.99132 7.977475 4.346715 0.8880713 13.964828
## 28.3 36.45846 9.705983 3.561876 0.7549162 16.677162
## 29.3 31.36729 9.026629 3.926023 0.7959465 15.416275
## 30.3 32.34048 7.943584 3.569362 0.8478048 13.917247# pairs(n0)
# plot(intervals(n0))
# Para conferir as predições.
pred <- with(tb_ref,
crossing(ue = levels(ue),
dia = seq(0, max(dia), by = 0.5)))
pred$y <- predict(n0, newdata = pred)
# Gráfico das curvas ajustadas sobre os valores observados.
ggplot(data = tb_ref,
mapping = aes(x = dia, y = nitro)) +
facet_wrap(facets = ~ue) +
geom_point(size = 1, pch = 3) +
geom_line(data = pred,
mapping = aes(x = dia, y = y),
color = "red") +
theme(strip.background = element_blank(),
strip.text = element_blank())
# Coefiente de determinação global.
cor(fitted(n0) + residuals(n0), fitted(n0))## [1] 0.9998445# Coeficiente de determinação por unidade experimental.
map_dbl(n0,
function(fit) {
cor(fitted(fit) + residuals(fit), fitted(fit))^2
}) %>%
sort()## 1.3 29.2 3.3 2.3 23.1 13.2 4.2
## 0.9987775 0.9990821 0.9991364 0.9992063 0.9992886 0.9993085 0.9993275
## 27.1 4.1 27.3 11.1 5.3 27.2 2.1
## 0.9993525 0.9993574 0.9993774 0.9993781 0.9993940 0.9994374 0.9994891
## 5.1 19.2 16.2 6.1 21.1 4.3 15.1
## 0.9994953 0.9995007 0.9995014 0.9995120 0.9995274 0.9995388 0.9995408
## 13.1 12.2 29.3 11.2 12.1 7.1 26.2
## 0.9995503 0.9995981 0.9996013 0.9996129 0.9996283 0.9996374 0.9996419
## 14.2 10.1 17.2 21.2 18.1 21.3 19.3
## 0.9996432 0.9996500 0.9996589 0.9996672 0.9996761 0.9996795 0.9996905
## 19.1 24.1 20.1 2.2 18.3 7.3 16.1
## 0.9996953 0.9996978 0.9996986 0.9997041 0.9997043 0.9997061 0.9997120
## 9.3 22.3 14.3 20.2 10.2 18.2 30.2
## 0.9997189 0.9997196 0.9997300 0.9997306 0.9997313 0.9997333 0.9997361
## 9.1 30.3 3.2 15.3 28.3 8.3 22.1
## 0.9997424 0.9997433 0.9997444 0.9997554 0.9997645 0.9997664 0.9997676
## 25.1 28.2 23.3 25.2 28.1 13.3 22.2
## 0.9997713 0.9997718 0.9997753 0.9997765 0.9997839 0.9997907 0.9997945
## 9.2 6.3 17.3 5.2 17.1 10.3 23.2
## 0.9997956 0.9997983 0.9998004 0.9998172 0.9998207 0.9998250 0.9998322
## 3.1 20.3 6.2 30.1 12.3 26.3 1.2
## 0.9998335 0.9998338 0.9998431 0.9998433 0.9998604 0.9998635 0.9998641
## 26.1 8.1 24.3 25.3 8.2 11.3 15.2
## 0.9998778 0.9998802 0.9998826 0.9998845 0.9998859 0.9998899 0.9998921
## 29.1 7.2 24.2 14.1 16.3 1.1
## 0.9998926 0.9998953 0.9999172 0.9999182 0.9999331 0.99993427 Análise das estimativas dos parâmetros como resposta
A análise considerando as estimativas dos parâmetros como variáveis respostas é uma alternativa.
#-----------------------------------------------------------------------
# Cria tabela com as estimativas dos parâmetros.
# tb_params_long <-
# summary(n0_fit)$coefficients %>%
# apply(MARGIN = 3,
# FUN = function(tb) {
# as.data.frame(tb) %>%
# rownames_to_column(var = "ue") %>%
# select(ue, Estimate, `Std. Error`) %>%
# rename("est" = "Estimate", "std" = "Std. Error")
# }) %>%
# bind_rows(.id = "param") %>%
# left_join(distinct(tb, ue, tipo_ureia, tipo_mo, dose_mo, bloco)) %>%
# mutate(Dose = factor(dose_mo)) %>%
# rename("Ureia" = "tipo_ureia", "Mo" = "tipo_mo", "dose" = dose_mo)
# head(tb_params_long)
# ggplot(data = tb_params_long,
# mapping = aes(x = Dose, y = std, color = Mo)) +
# facet_wrap(facets = Ureia ~ param, scale = "free") +
# geom_point()
# Cria tabela com as estimativas dos parâmetros e as condições
# experimetais.
tb_params <-
coef(n0) %>%
as.data.frame() %>%
rownames_to_column(var = "ue") %>%
left_join(distinct(tb, ue, tipo_ureia, tipo_mo, dose_mo, bloco)) %>%
mutate(Dose = factor(dose_mo)) %>%
rename("Ureia" = "tipo_ureia", "Mo" = "tipo_mo", "dose" = dose_mo)## Joining, by = "ue"#-----------------------------------------------------------------------
# Calcula a taxa máxima no ponto de inflexão.
# 1ª derivada da função.
D(expression(tha/(1 + (thv/x)^thk)), "x")## tha * ((thv/x)^(thk - 1) * (thk * (thv/x^2)))/(1 + (thv/x)^thk)^2# Função que determina o ponto de inflexão e a taxa nele.
calc_params <- function(params) {
obj_fun <- function(x, params) {
with(params, {
tha * ((thv/x)^(thk - 1) *
(thk * (thv/x^2)))/(1 + (thv/x)^thk)^2
})
}
opt <- optimise(f = obj_fun,
params,
interval = c(0, 25),
maximum = TRUE)
return(data.frame(thi = opt$maximum, thr = opt$objective))
}
calc_params(start)## thi thr
## 1 7.143304 3.173134# Calcula e adiciona as estimativas das funções de parâmetros.
tb_params <- tb_params %>%
rowwise() %>%
do({
parlist <- as.list(.[names(coef(n0))])
cbind(., calc_params(parlist))
}) %>%
ungroup()
lattice::splom(tb_params[, grep(x = names(tb_params), "^th")],
# group = tb_params$Ureia
# group = tb_params$Mo
# group = tb_params$Dose
)
7.1 Conteúdo de N de fácil liberação
#-----------------------------------------------------------------------
# theta_a: assíntota do conteúdo de pronta liberação.
# Modelo maximal onde todos os termos são qualitativos.
m0 <- lm(tha ~ bloco + Ureia * Mo * Dose,
data = tb_params)
# m0 <- lm(est ~ bloco + Ureia * Mo * Dose,
# weights = 1/(std^2),
# data = filter(tb_params_long, param == "A"))
# Modelo em que dose entre com polinômio de 2º grau.
m1 <- update(m0,
formula = . ~ bloco +
Ureia * Mo * poly(dose, degree = 2))
# Teste da falta de ajuste ou significância dos termos abandonados.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: tha ~ bloco + Ureia * Mo * Dose
## Model 2: tha ~ bloco + Ureia + Mo + poly(dose, degree = 2) + Ureia:Mo +
## Ureia:poly(dose, degree = 2) + Mo:poly(dose, degree = 2) +
## Ureia:Mo:poly(dose, degree = 2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 1189.9
## 2 70 1734.3 -12 -544.35 2.2111 0.02277 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Diagnóstico dos resíduos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: tha
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 2 11.13 5.567 0.2714 0.7632937
## Ureia 2 148.15 74.077 3.6107 0.0332730 *
## Mo 1 45.03 45.026 2.1947 0.1438973
## Dose 4 530.86 132.715 6.4689 0.0002247 ***
## Ureia:Mo 2 4.58 2.288 0.1115 0.8946533
## Ureia:Dose 8 429.39 53.674 2.6162 0.0162039 *
## Mo:Dose 4 77.75 19.438 0.9475 0.4432450
## Ureia:Mo:Dose 8 112.72 14.090 0.6868 0.7013767
## Residuals 58 1189.91 20.516
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Comparações múltiplas.
tb_emm <- emmeans(m0, specs = ~Dose | Ureia) %>%
multcomp::cld(Letters = letters)## NOTE: Results may be misleading due to involvement in interactionstb_emm## Ureia = incorporado:
## Dose emmean SE df lower.CL upper.CL .group
## 0 21.3 1.85 58 17.6 25.0 a
## 150 31.5 1.85 58 27.8 35.2 b
## 600 31.6 1.85 58 27.9 35.3 b
## 1200 34.2 1.85 58 30.5 37.9 b
## 300 35.9 1.85 58 32.2 39.6 b
##
## Ureia = revestido:
## Dose emmean SE df lower.CL upper.CL .group
## 0 28.2 1.85 58 24.5 31.9 a
## 300 31.0 1.85 58 27.3 34.7 a
## 600 32.6 1.85 58 28.9 36.3 a
## 1200 34.2 1.85 58 30.5 37.9 a
## 150 34.2 1.85 58 30.5 37.9 a
##
## Ureia = ureia:
## Dose emmean SE df lower.CL upper.CL .group
## 0 33.3 1.85 58 29.6 37.0 a
## 600 33.3 1.85 58 29.6 37.0 a
## 1200 33.6 1.85 58 29.9 37.3 a
## 300 34.0 1.85 58 30.3 37.7 a
## 150 35.8 1.85 58 32.1 39.6 a
##
## Results are averaged over the levels of: bloco, Mo
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05# Prepara a tabela.
tb_emm <- tb_emm %>%
as.data.frame() %>%
rename("lwr" = 6, "upr" = 7, "cld" = ".group") %>%
mutate(cld = str_trim(cld))
# Gráfico com os resultados.
ggplot(data = tb_emm,
mapping = aes(x = Dose, y = emmean, ymin = lwr, ymax = upr)) +
facet_wrap(facets = ~Ureia) +
geom_errorbar(width = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.1f%s", emmean, cld)),
vjust = -0.25) +
labs(x = "Dose de molibdênio",
y = "Conteúdo de N de pronta liberação")
7.2 Tempo de meia-vida ou EC50
#-----------------------------------------------------------------------
# theta_v: tempo de meia vida ou EC50.
# Modelo maximal onde todos os termos são qualitativos.
m0 <- lm(thv ~ bloco + Ureia * Mo * Dose,
data = tb_params)
# m0 <- lm(est ~ bloco + Ureia * Mo * Dose,
# weights = 1/(std^2),
# data = filter(tb_params_long, param == "R"))
# Modelo em que dose entre com polinômio de 2º grau.
m1 <- update(m0,
formula = . ~ bloco +
Ureia * Mo * poly(dose, degree = 2))
# Teste da falta de ajuste ou significância dos termos abandonados.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: thv ~ bloco + Ureia * Mo * Dose
## Model 2: thv ~ bloco + Ureia + Mo + poly(dose, degree = 2) + Ureia:Mo +
## Ureia:poly(dose, degree = 2) + Mo:poly(dose, degree = 2) +
## Ureia:Mo:poly(dose, degree = 2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 78.588
## 2 70 124.606 -12 -46.017 2.8301 0.004067 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Diagnóstico dos resíduos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# MASS::boxcox(m0)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: thv
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 2 4.776 2.3881 1.7625 0.18069
## Ureia 2 5.098 2.5489 1.8812 0.16159
## Mo 1 3.916 3.9155 2.8898 0.09450 .
## Dose 4 47.549 11.8873 8.7731 1.311e-05 ***
## Ureia:Mo 2 3.601 1.8006 1.3289 0.27271
## Ureia:Dose 8 41.182 5.1477 3.7991 0.00121 **
## Mo:Dose 4 6.683 1.6708 1.2331 0.30690
## Ureia:Mo:Dose 8 10.214 1.2768 0.9423 0.48945
## Residuals 58 78.588 1.3550
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Comparações múltiplas.
tb_emm <- emmeans(m0, specs = ~Dose | Ureia) %>%
multcomp::cld(Letters = letters)## NOTE: Results may be misleading due to involvement in interactionstb_emm## Ureia = incorporado:
## Dose emmean SE df lower.CL upper.CL .group
## 150 8.42 0.475 58 7.47 9.37 a
## 1200 8.59 0.475 58 7.63 9.54 a
## 600 8.63 0.475 58 7.68 9.59 a
## 300 8.93 0.475 58 7.98 9.88 a
## 0 12.50 0.475 58 11.55 13.45 b
##
## Ureia = revestido:
## Dose emmean SE df lower.CL upper.CL .group
## 1200 8.59 0.475 58 7.64 9.54 a
## 600 8.62 0.475 58 7.67 9.57 a
## 150 8.91 0.475 58 7.96 9.86 a
## 300 8.97 0.475 58 8.02 9.92 a
## 0 10.10 0.475 58 9.15 11.05 a
##
## Ureia = ureia:
## Dose emmean SE df lower.CL upper.CL .group
## 600 8.20 0.475 58 7.25 9.15 a
## 1200 8.45 0.475 58 7.50 9.41 a
## 0 8.81 0.475 58 7.86 9.76 a
## 150 9.13 0.475 58 8.18 10.08 a
## 300 9.60 0.475 58 8.65 10.56 a
##
## Results are averaged over the levels of: bloco, Mo
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05# Prepara a tabela.
tb_emm <- tb_emm %>%
as.data.frame() %>%
rename("lwr" = 6, "upr" = 7, "cld" = ".group") %>%
mutate(cld = str_trim(cld))
# Gráfico com os resultados.
ggplot(data = tb_emm,
mapping = aes(x = Dose, y = emmean, ymin = lwr, ymax = upr)) +
facet_wrap(facets = ~Ureia) +
geom_errorbar(width = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.2f%s", emmean, cld)),
vjust = -0.25) +
labs(x = "Dose de molibdênio",
y = "Período para liberação de 50% do N de pronta liberação")
7.3 Taxa máxima de acúmulo de N
#-----------------------------------------------------------------------
# theta_r: taxa na inflexão ou taxa máxima da fase pronta.
# Modelo maximal onde todos os termos são qualitativos.
m0 <- lm(thr ~ bloco + Ureia * Mo * Dose,
data = tb_params)
# m0 <- lm(est ~ bloco + Ureia * Mo * Dose,
# weights = 1/(std^2),
# data = filter(tb_params_long, param == "R"))
# Modelo em que dose entre com polinômio de 2º grau.
m1 <- update(m0,
formula = . ~ bloco +
Ureia * Mo * poly(dose, degree = 2))
# Teste da falta de ajuste ou significância dos termos abandonados.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: thr ~ bloco + Ureia * Mo * Dose
## Model 2: thr ~ bloco + Ureia + Mo + poly(dose, degree = 2) + Ureia:Mo +
## Ureia:poly(dose, degree = 2) + Mo:poly(dose, degree = 2) +
## Ureia:Mo:poly(dose, degree = 2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 21.101
## 2 70 37.222 -12 -16.121 3.6927 0.0003795 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Diagnóstico dos resíduos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# MASS::boxcox(m0)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: thr
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 2 0.9160 0.45802 1.2590 0.291594
## Ureia 2 2.1023 1.05117 2.8894 0.063651 .
## Mo 1 0.1695 0.16946 0.4658 0.497646
## Dose 4 11.5028 2.87570 7.9045 3.719e-05 ***
## Ureia:Mo 2 0.3196 0.15981 0.4393 0.646633
## Ureia:Dose 8 10.6273 1.32841 3.6514 0.001666 **
## Mo:Dose 4 3.5910 0.89776 2.4677 0.054713 .
## Ureia:Mo:Dose 8 5.4082 0.67602 1.8582 0.084510 .
## Residuals 58 21.1008 0.36381
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Comparações múltiplas.
tb_emm <- emmeans(m0, specs = ~Dose | Ureia) %>%
multcomp::cld(Letters = letters)## NOTE: Results may be misleading due to involvement in interactionstb_emm## Ureia = incorporado:
## Dose emmean SE df lower.CL upper.CL .group
## 0 1.63 0.246 58 1.14 2.12 a
## 600 3.31 0.246 58 2.82 3.81 b
## 300 3.56 0.246 58 3.06 4.05 b
## 150 3.56 0.246 58 3.07 4.05 b
## 1200 3.64 0.246 58 3.15 4.13 b
##
## Ureia = revestido:
## Dose emmean SE df lower.CL upper.CL .group
## 0 2.69 0.246 58 2.19 3.18 a
## 600 3.19 0.246 58 2.70 3.68 ab
## 300 3.30 0.246 58 2.81 3.80 ab
## 150 3.40 0.246 58 2.91 3.89 ab
## 1200 3.68 0.246 58 3.19 4.17 b
##
## Ureia = ureia:
## Dose emmean SE df lower.CL upper.CL .group
## 300 3.10 0.246 58 2.60 3.59 a
## 150 3.50 0.246 58 3.01 3.99 a
## 0 3.58 0.246 58 3.09 4.07 a
## 600 3.61 0.246 58 3.12 4.10 a
## 1200 3.75 0.246 58 3.26 4.24 a
##
## Results are averaged over the levels of: bloco, Mo
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 5 estimates
## significance level used: alpha = 0.05# Prepara a tabela.
tb_emm <- tb_emm %>%
as.data.frame() %>%
rename("lwr" = 6, "upr" = 7, "cld" = ".group") %>%
mutate(cld = str_trim(cld))
# Gráfico com os resultados.
ggplot(data = tb_emm,
mapping = aes(x = Dose, y = emmean, ymin = lwr, ymax = upr)) +
facet_wrap(facets = ~Ureia) +
geom_errorbar(width = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.2f%s", emmean, cld)),
vjust = -0.25) +
labs(x = "Dose de molibdênio",
y = "Taxa máxima de acúmulo de N de pronta liberação")
7.4 Período para início da fase lenta
#-----------------------------------------------------------------------
# theta_b: período para início da fase de lenta liberação.
# Modelo maximal onde todos os termos são qualitativos.
m0 <- lm(thb ~ bloco + Ureia * Mo * Dose,
data = tb_params)
# m0 <- lm(est ~ bloco + Ureia * Mo * Dose,
# weights = 1/(std^2),
# data = filter(tb_params_long, param == "K"))
# Modelo em que dose entre com polinômio de 2º grau.
m1 <- update(m0,
formula = . ~ bloco +
Ureia * Mo * poly(dose, degree = 2))
# Teste da falta de ajuste ou significância dos termos abandonados.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: thb ~ bloco + Ureia * Mo * Dose
## Model 2: thb ~ bloco + Ureia + Mo + poly(dose, degree = 2) + Ureia:Mo +
## Ureia:poly(dose, degree = 2) + Mo:poly(dose, degree = 2) +
## Ureia:Mo:poly(dose, degree = 2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 458.26
## 2 70 554.67 -12 -96.408 1.0168 0.446# Diagnóstico dos resíduos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: thb
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 2 9.92 4.958 0.6276 0.53748
## Ureia 2 48.87 24.437 3.0929 0.05293 .
## Mo 1 35.41 35.406 4.4812 0.03857 *
## Dose 4 65.60 16.401 2.0758 0.09565 .
## Ureia:Mo 2 30.23 15.113 1.9128 0.15687
## Ureia:Dose 8 59.48 7.435 0.9410 0.49044
## Mo:Dose 4 9.89 2.473 0.3130 0.86819
## Ureia:Mo:Dose 8 32.07 4.009 0.5074 0.84593
## Residuals 58 458.26 7.901
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Comparações múltiplas.
tb_emm <- emmeans(m0, specs = ~Mo) %>%
multcomp::cld(Letters = letters)## NOTE: Results may be misleading due to involvement in interactionstb_emm## Mo emmean SE df lower.CL upper.CL .group
## tmo 13.8 0.419 58 13.0 14.6 a
## ma 15.0 0.419 58 14.2 15.9 b
##
## Results are averaged over the levels of: bloco, Ureia, Dose
## Confidence level used: 0.95
## significance level used: alpha = 0.05# Prepara a tabela.
tb_emm <- tb_emm %>%
as.data.frame() %>%
rename("lwr" = 5, "upr" = 6, "cld" = ".group") %>%
mutate(cld = str_trim(cld))
# Gráfico com os resultados.
ggplot(data = tb_emm,
mapping = aes(x = Mo, y = emmean, ymin = lwr, ymax = upr)) +
geom_errorbar(width = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.2f%s", emmean, cld)),
vjust = -0.25) +
labs(x = "Fonte de molibdênio",
y = "Período para início da fase de lenta liberação")
7.5 Taxa de acúmulo de lenta liberação
#-----------------------------------------------------------------------
# theta_x: taxa da fase de lenta liberação.
# Modelo maximal onde todos os termos são qualitativos.
m0 <- lm(thx ~ bloco + Ureia * Mo * Dose,
data = tb_params)
# m0 <- lm(est ~ bloco + Ureia * Mo * Dose,
# weights = 1/(std^2),
# data = filter(tb_params_long, param == "B"))
# Modelo em que dose entre com polinômio de 2º grau.
m1 <- update(m0,
formula = . ~ bloco +
Ureia * Mo * poly(dose, degree = 2))
# Teste da falta de ajuste ou significância dos termos abandonados.
anova(m0, m1)## Analysis of Variance Table
##
## Model 1: thx ~ bloco + Ureia * Mo * Dose
## Model 2: thx ~ bloco + Ureia + Mo + poly(dose, degree = 2) + Ureia:Mo +
## Ureia:poly(dose, degree = 2) + Mo:poly(dose, degree = 2) +
## Ureia:Mo:poly(dose, degree = 2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 58 1.7290
## 2 70 2.0331 -12 -0.30412 0.8501 0.6# Diagnóstico dos resíduos.
par(mfrow = c(2, 2))
plot(m0)
layout(1)
# Quadro de análise de variância.
anova(m0)## Analysis of Variance Table
##
## Response: thx
## Df Sum Sq Mean Sq F value Pr(>F)
## bloco 2 0.02159 0.010793 0.3620 0.69781
## Ureia 2 0.00948 0.004742 0.1591 0.85330
## Mo 1 0.01988 0.019884 0.6670 0.41744
## Dose 4 0.00431 0.001076 0.0361 0.99744
## Ureia:Mo 2 0.17840 0.089199 2.9922 0.05798 .
## Ureia:Dose 8 0.25099 0.031374 1.0524 0.40866
## Mo:Dose 4 0.01502 0.003756 0.1260 0.97249
## Ureia:Mo:Dose 8 0.12856 0.016071 0.5391 0.82221
## Residuals 58 1.72902 0.029811
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Comparações múltiplas.
tb_emm <- emmeans(m0, specs = ~Mo | Ureia) %>%
multcomp::cld(Letters = letters)## NOTE: Results may be misleading due to involvement in interactionstb_emm## Ureia = incorporado:
## Mo emmean SE df lower.CL upper.CL .group
## tmo 0.734 0.0446 58 0.645 0.823 a
## ma 0.814 0.0446 58 0.725 0.903 a
##
## Ureia = revestido:
## Mo emmean SE df lower.CL upper.CL .group
## ma 0.689 0.0446 58 0.600 0.778 a
## tmo 0.827 0.0446 58 0.738 0.916 b
##
## Ureia = ureia:
## Mo emmean SE df lower.CL upper.CL .group
## ma 0.734 0.0446 58 0.644 0.823 a
## tmo 0.765 0.0446 58 0.675 0.854 a
##
## Results are averaged over the levels of: bloco, Dose
## Confidence level used: 0.95
## significance level used: alpha = 0.05# Prepara a tabela.
tb_emm <- tb_emm %>%
as.data.frame() %>%
rename("lwr" = 6, "upr" = 7, "cld" = ".group") %>%
mutate(cld = str_trim(cld))
# Gráfico com os resultados.
ggplot(data = tb_emm,
mapping = aes(x = Mo, y = emmean, ymin = lwr, ymax = upr)) +
facet_wrap(facets = ~Ureia) +
geom_errorbar(width = 0.1) +
geom_point() +
geom_label(mapping = aes(label = sprintf("%0.2f%s", emmean, cld)),
vjust = -0.25) +
labs(x = "Fonte de molibdênio",
y = "Taxa de acúmulo de N de difícil liberação")