Ajuste de curvas de progresso de doença em folhas
Henrique da Silva Silveira Duarte e Walmes Marques Zeviani
2019-11-01
Leitura e inspeção dos dados
rm(list = objects())
library(tidyverse)
library(nlme)
library(multcomp)
## Importa os dados.
tb <- read_tsv("pustulas.txt", locale = locale(decimal_mark = ","))
attr(tb, "spec") <- NULL
# Renomeia.
names(tb) <- c("estado", "rept", "folha", "dia", "pustulas")
# Som o dia em que apareceu o primeiro sintoma em cada estado.
tb <- tb %>%
mutate(dia_u = dia,
dia = case_when(estado == "SP" ~ dia_u + 9,
estado == "SC" ~ dia_u + 7,
estado == "PR" ~ dia_u + 10))
# Filtra, converte e cria variáveis.
tb <- filter(tb,
estado %in% c("SP", "PR", "SC"),
dia_u > 0) %>%
mutate(estado = factor(estado),
rept = factor(rept),
folha = factor(folha),
rept_folha = interaction(rept, folha, drop = TRUE),
ue = interaction(estado, rept, folha, drop = TRUE))
str(tb)## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 288 obs. of 8 variables:
## $ estado : Factor w/ 3 levels "PR","SC","SP": 3 3 3 3 3 3 3 3 3 3 ...
## $ rept : Factor w/ 4 levels "1","2","3","4": 1 1 2 2 3 3 4 4 1 1 ...
## $ folha : Factor w/ 2 levels "1","2": 1 2 1 2 1 2 1 2 1 2 ...
## $ dia : num 10 10 10 10 10 10 10 10 11 11 ...
## $ pustulas : num 18 2 1 16 28 4 11 2 29 4 ...
## $ dia_u : num 1 1 1 1 1 1 1 1 2 2 ...
## $ rept_folha: Factor w/ 8 levels "1.1","2.1","3.1",..: 1 5 2 6 3 7 4 8 1 5 ...
## $ ue : Factor w/ 24 levels "PR.1.1","SC.1.1",..: 3 15 6 18 9 21 12 24 3 15 ...# Tabela de frequência.
xtabs(~estado + rept, data = tb)## rept
## estado 1 2 3 4
## PR 22 22 22 22
## SC 26 26 26 26
## SP 24 24 24 24# Exploratória.
ggplot(data = tb,
mapping = aes(x = dia,
y = pustulas,
group = rept_folha,
color = rept)) +
facet_wrap(facets = ~estado, nrow = 1) +
expand_limits(x = c(0, NA)) +
geom_point() +
geom_line()
ggplot(data = tb,
mapping = aes(x = dia,
y = pustulas,
group = folha,
color = folha)) +
facet_grid(facets = estado ~ rept) +
# scale_y_log10() +
expand_limits(x = c(0, NA)) +
geom_point() +
geom_line()
# ggplot(data = tb,
# mapping = aes(x = dia,
# y = pustulas)) +
# facet_wrap(facets = ~ue) +
# geom_point() +
# geom_line()
tb <- tb %>%
group_by(estado, rept, dia) %>%
summarise(pustulas = mean(pustulas)) %>%
ungroup() %>%
mutate(ue = interaction(estado, rept, drop = TRUE))
#Ajuste das curvas para o progresso da doença
Foi usado o modelo logÃstico na parametrização centro-escala \[ f(x, \theta_a, \theta_c, \theta_s) = \dfrac{\theta_a}{1 + \exp\left\{-\dfrac{x - \theta_c}{\theta_s}\right\}} \] em que:
- \(\theta_a\) é a assÃntota do modelo e representa o limite da função quando o tempo vai para o infinito (\(\lim_{x \to \infty} f(x) = \theta_a\)).
- \(\theta_c\) é o valor no tempo correspondente ao ponto de inflexão, ponto em que a taxa de crescimento é máxima e no qual houve metade do crescimento máximo \(\theta_a\), ou seja, \(f(x = \theta_c) = \theta_{a}/2\).
- \(\theta_s\) é o parâmetro de escala relativa. Ele se relaciona de forma diretamente inversa com a taxa no ponto de inflexão (\(\theta_r\)) por \(\theta_r = \theta_a/(4 \theta_s)\). Dessa forma, o significado é \(\theta_s = \theta_a/(4 \theta_s)\), portanto é 1/4 da taxa relativa de crescimento.
Para haver uma interpração direta do modelo em termos dos parâmetros \(\theta_a\), \(\theta_c\) e \(\theta_r\), que tem maior apelo biológico, fez-se a reparametrização do modelo logÃstico substituindo-se \(\theta_s\) por \(\theta_r\),
\[ f(x, \theta_a, \theta_c, \theta_r) = \dfrac{\theta_a}{1 + \exp\left\{-4\theta_r\dfrac{(x - \theta_c)}{\theta_a}\right\}}. \]
# Ajuste de uma curva a todos os dados (passo 1 para ganhar experiência
# com o modelo até chegar ao modelo adequado).
plot(pustulas ~ dia, data = tb)
start <- list(A = 150, C = 16, R = 10)
m0 <- nls(pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)),
start = start,
data = tb)
coef(m0)## A C R
## 138.22495 14.64895 15.00900# Cria estrutura de dados agrupados.
db <- groupedData(pustulas ~ dia | ue,
data = tb,
order.groups = FALSE)
str(db)## Classes 'nfnGroupedData', 'nfGroupedData', 'groupedData', 'tbl_df', 'tbl' and 'data.frame': 144 obs. of 5 variables:
## $ estado : Factor w/ 3 levels "PR","SC","SP": 1 1 1 1 1 1 1 1 1 1 ...
## $ rept : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
## $ dia : num 11 12 13 14 15 16 17 18 19 20 ...
## $ pustulas: num 0.5 3 7 12.5 23 30 38.5 48.5 54.5 57.5 ...
## $ ue : Factor w/ 12 levels "PR.1","SC.1",..: 1 1 1 1 1 1 1 1 1 1 ...
## - attr(*, "formula")=Class 'formula' language pustulas ~ dia | ue
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## - attr(*, "FUN")=function (x)
## - attr(*, "order.groups")= logi FALSE# Ajuste de uma curva para cada estado que ignora o efeito de
# repetição:folha (passo 2 para ganhar experiência com o modelo até
# chegar ao modelo adequado).
m0 <- nlsList(pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)) | estado,
start = start,
data = db)
coef(m0)## A C R
## PR 67.04408 16.16532 8.717485
## SC 214.24142 13.82643 23.932080
## SP 204.01890 17.55130 18.241930plot(augPred(m0))
# Usa as estimativas para passar como valores iniciais no modelo final.
coefs <- as.matrix(coef(m0))
start <- rbind(coefs[1, , drop = FALSE],
sweep(x = coefs[-1, , drop = FALSE],
MARGIN = 2,
STATS = coefs[1, , drop = FALSE],
FUN = "-"))
c(start)## [1] 67.044076 147.197348 136.974821 16.165320 -2.338894 1.385975
## [7] 8.717485 15.214595 9.524445# Modelo com efeito fixo de estado nos parâmetros da curva e efeito
# aleatório da unidade experimental nos mesmos parâmetros.
mm0 <- nlme(model = pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)),
data = db,
fixed = A + C + R ~ estado,
random = A + C + R ~ 1,
start = c(start))
# Resumo do ajuste.
summary(mm0)## Nonlinear mixed-effects model fit by maximum likelihood
## Model: pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A))
## Data: db
## AIC BIC logLik
## 927.9684 975.4854 -447.9842
##
## Random effects:
## Formula: list(A ~ 1, C ~ 1, R ~ 1)
## Level: ue
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## A.(Intercept) 48.2816632 A.(In) C.(In)
## C.(Intercept) 0.4000299 0.159
## R.(Intercept) 4.2469459 0.931 0.163
## Residual 3.6974245
##
## Fixed effects: A + C + R ~ estado
## Value Std.Error DF t-value p-value
## A.(Intercept) 67.30344 25.09122 124 2.68235 0.0083
## A.estadoSC 147.33060 35.44537 124 4.15655 0.0001
## A.estadoSP 133.01875 36.03918 124 3.69095 0.0003
## C.(Intercept) 16.16015 0.31855 124 50.73044 0.0000
## C.estadoSC -2.35220 0.38781 124 -6.06527 0.0000
## C.estadoSP 1.27573 0.44181 124 2.88753 0.0046
## R.(Intercept) 8.78100 2.24639 124 3.90893 0.0002
## R.estadoSC 15.52146 3.16666 124 4.90152 0.0000
## R.estadoSP 9.60549 3.15438 124 3.04512 0.0028
## Correlation:
## A.(In) A.stSC A.stSP C.(In) C.stSC C.stSP R.(In) R.stSC
## A.estadoSC -0.708
## A.estadoSP -0.696 0.493
## C.(Intercept) 0.179 -0.127 -0.125
## C.estadoSC -0.147 0.175 0.103 -0.821
## C.estadoSP -0.129 0.092 0.238 -0.721 0.592
## R.(Intercept) 0.892 -0.632 -0.621 0.034 -0.028 -0.025
## R.estadoSC -0.633 0.898 0.441 -0.024 0.070 0.018 -0.709
## R.estadoSP -0.635 0.450 0.892 -0.024 0.020 0.073 -0.712 0.505
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.67482017 -0.57379992 -0.06831882 0.49403889 3.09335607
##
## Number of Observations: 144
## Number of Groups: 12# Quadro de testes de hipótese para o efeito de estado nos parâmetros da
# curva.
anova(mm0, type = "marginal")## numDF denDF F-value p-value
## A.(Intercept) 1 124 7.1950 0.0083
## A.estado 2 124 10.4200 0.0001
## C.(Intercept) 1 124 2573.5772 <.0001
## C.estado 2 124 50.7276 <.0001
## R.(Intercept) 1 124 15.2797 0.0002
## R.estado 2 124 12.2297 <.0001# Curvas ajustadas para cada unidade experimental.
plot(augPred(mm0))
# Gráfico de pares dos valores preditos dos efeitos aleatórios.
pairs(mm0, as.matrix = TRUE)
# Ajusta o modelo com estimativas por ponto experimental.
mm1 <- nlme(model = pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)),
data = db,
fixed = A + C + R ~ 0 + estado,
random = A + C + R ~ 1,
start = unlist(coef(m0)))
fixef(mm1)## A.estadoPR A.estadoSC A.estadoSP C.estadoPR C.estadoSC C.estadoSP
## 67.303442 214.634041 200.322191 16.160151 13.807952 17.435876
## R.estadoPR R.estadoSC R.estadoSP
## 8.780996 24.302453 18.386491# Verifica se as log-verossimilhanças são iguais.
logLik(mm0)## 'log Lik.' -447.9842 (df=16)logLik(mm1)## 'log Lik.' -447.9842 (df=16)# Extrai os intervalos de confiança para as estimativas de efeito fixo.
ci <- intervals(mm1, which = "fixed")
ci <- as.data.frame(ci$fixed) %>%
rownames_to_column(var = "caso") %>%
separate("caso", c("param", "estado"), sep = "\\.estado") %>%
rename("fit" = "est.")
ci## param estado lower fit upper
## 1 A PR 19.217878 67.303442 115.38901
## 2 A SC 166.654181 214.634041 262.61390
## 3 A SP 150.744308 200.322191 249.90007
## 4 C PR 15.549674 16.160151 16.77063
## 5 C SC 13.384052 13.807952 14.23185
## 6 C SP 16.849190 17.435876 18.02256
## 7 R PR 4.475939 8.780996 13.08605
## 8 R SC 20.025138 24.302453 28.57977
## 9 R SP 14.142625 18.386491 22.63036#Testes de hipótese para os parâmetros da curva
# fixef(mm1)
# vcov(mm1)
# Define os contrastes.
K <- rbind(c(1, -1, 0), # H_0: PR - SC.
c(1, -1, -1), # H_0: PR - SP.
c(0, 1, -1)) # H_0: SC - SP.
rownames(K) <- c("PR-SC", "PR-SP", "SC-SP")
# Matrizes para contrastes par a par.
K_A <- cbind(K, diag(3) * 0, diag(3) * 0)
K_C <- cbind(diag(3) * 0, K, diag(3) * 0)
K_R <- cbind(diag(3) * 0, diag(3) * 0, K)
# Testes sobre o parâmetro de assÃntota A.
summary(glht(model = mm1, linfct = K_A))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: nlme.formula(model = pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)),
## data = db, fixed = A + C + R ~ 0 + estado, random = A + C +
## R ~ 1, start = unlist(coef(m0)))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## PR-SC == 0 -147.33 34.32 -4.293 <1e-04 ***
## PR-SP == 0 -347.65 42.49 -8.182 <1e-04 ***
## SC-SP == 0 14.31 34.86 0.411 0.926
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# Testes sobre o parâmetro de centro C.
summary(glht(model = mm1, linfct = K_C))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: nlme.formula(model = pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)),
## data = db, fixed = A + C + R ~ 0 + estado, random = A + C +
## R ~ 1, start = unlist(coef(m0)))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## PR-SC == 0 2.3522 0.3755 6.264 7.32e-10 ***
## PR-SP == 0 -15.0837 0.4784 -31.530 < 1e-10 ***
## SC-SP == 0 -3.6279 0.3657 -9.921 < 1e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# Testes sobre o parâmetro de taxa R.
summary(glht(model = mm1, linfct = K_R))##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: nlme.formula(model = pustulas ~ A/(1 + exp(-4 * R * (dia - C)/A)),
## data = db, fixed = A + C + R ~ 0 + estado, random = A + C +
## R ~ 1, start = unlist(coef(m0)))
##
## Linear Hypotheses:
## Estimate Std. Error z value Pr(>|z|)
## PR-SC == 0 -15.521 3.066 -5.062 <0.001 ***
## PR-SP == 0 -33.908 3.741 -9.063 <0.001 ***
## SC-SP == 0 5.916 3.044 1.943 0.121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)# Adiciona as letras manualmente com base nos p-valores.
ci$let <- c("b", "a", "a",
"b", "c", "a",
"b", "a", "a")
# Gráfico com segmentos.
gg_ci <-
ggplot(data = ci,
mapping = aes(x = estado, y = fit)) +
facet_wrap(facets = ~param,
scale = "free_y",
labeller = labeller(param = c(A = "Asymptote (A)",
C = "Center (C)",
R = "Rate (R)"))) +
geom_point() +
geom_errorbar(mapping = aes(ymin = lower, ymax = upper),
width = 0.1) +
# geom_text(mapping = aes(label = sprintf("%s", let)),
# hjust = 0,
# nudge_x = 0.075) +
geom_text(mapping = aes(label = sprintf("%0.2f %s",
fit,
let)),
hjust = 0.5,
vjust = -1,
angle = 90,
nudge_x = 0.075) +
labs(x = "Origin",
y = "Estimate with confidence interval (95%)")
gg_ci
#Curvas ajustadas com bandas de confiança
# help(SSlogis, h = "html")
model <- ~A/(1 + exp(-4 * R * (dia - C)/A))
# Para retornar a matriz gradiente..
partials <- deriv(model,
namevec = c("A", "C", "R"),
function.arg = function(dia, A, C, R) {
NULL
})
# Estimativas pontuais para cada nÃvel do efeito fixo.
coefs <- matrix(fixef(mm1), nrow = 3)
# Matriz com os parâmetros estimados por solo.
params <- matrix(coefs,
ncol = 3,
dimnames = list(levels(db$estado),
c("A", "C", "R")))
# Cria grid de pontos para predição.
dia_seq <- seq(0, max(tb$dia), by = 0.25)
grid <- with(db,
expand.grid(dia = dia_seq,
estado = levels(estado)))
grid$fit <- predict(mm1, newdata = grid, level = 0)
# Lista de matrizes gradiente.
grad <- sapply(levels(db$estado),
simplify = FALSE,
FUN = function(nome) {
u <- do.call(partials,
args = c(list(dia = dia_seq),
as.list(params[nome, ])))
u <- attr(u, "gradient")
colnames(u) <- sprintf("%s.estado%s",
colnames(u),
nome)
return(u)
})
# Cria matriz bloco diagonal e atribui nomes.
X <- as.matrix(Matrix::bdiag(grad))
colnames(X) <- c(sapply(grad, colnames))
# Troca colunas de lugar para corresponder com `vcov()`.
X <- X[, colnames(vcov(mm1))]
head(X)## A.estadoPR A.estadoSC A.estadoSP C.estadoPR C.estadoSC C.estadoSP
## [1,] 0.002050410 0 0 -0.007634102 0 0
## [2,] 0.002303719 0 0 -0.008697482 0 0
## [3,] 0.002587792 0 0 -0.009908899 0 0
## [4,] 0.002906281 0 0 -0.011288939 0 0
## [5,] 0.003263254 0 0 -0.012861039 0 0
## [6,] 0.003663246 0 0 -0.014651888 0 0
## R.estadoPR R.estadoSC R.estadoSP
## [1,] -0.01404946 0 0
## [2,] -0.01575883 0 0
## [3,] -0.01767167 0 0
## [4,] -0.01981145 0 0
## [5,] -0.02220423 0 0
## [6,] -0.02487894 0 0# Cria erro padrão do valor predito e calcula o intervalo de confiança.
qtl <- qt(0.975, df = nrow(db) - length(fixef(mm1)))
grid$se <- sqrt(diag(X %*% vcov(mm1) %*% t(X)))
grid <- grid %>%
mutate(lwr = fit - qtl * se,
upr = fit + qtl * se)
# Curvas ajustadas com bandas de confiança para a média.
ggplot(data = grid,
mapping = aes(x = dia, y = fit)) +
facet_wrap(facets = ~estado, nrow = 1) +
geom_ribbon(data = grid,
mapping = aes(ymin = lwr, ymax = upr),
fill = "darkorange",
alpha = 0.4) +
geom_line()
# Tabela com texto para ser exibido no gráfico.
params <- params %>%
as.data.frame() %>%
rownames_to_column(var = "estado")
fmt <- "Asymp: %0.2f\nCenter: %0.2f\nRate: %0.2f"
params$label <- with(params, sprintf(fmt, A, C, R))
# Curvas ajustadas com bandas de confiança para a média e valores
# observados.
gg_cb <-
ggplot(data = tb,
mapping = aes(x = dia,
y = pustulas,
# group = rept_folha)) +
group = ue)) +
facet_wrap(facets = ~estado,
nrow = 1) +
geom_point(color = "gray50",
pch = 1) +
geom_line(color = "gray50",
size = 0.2) +
geom_ribbon(data = grid,
inherit.aes = FALSE,
mapping = aes(x = dia, ymin = lwr, ymax = upr),
fill = "darkorange",
alpha = 0.4) +
geom_line(data = grid,
inherit.aes = FALSE,
mapping = aes(x = dia, y = fit)) +
geom_text(data = params,
inherit.aes = FALSE,
mapping = aes(x = 0.5, y = 350, label = label),
parse = FALSE,
size = 3.5,
hjust = 0,
vjust = 1) +
labs(x = "Days after inoculation",
y = expression("Number of pustules" ~ cm^{-2}))
gg_cb
#gridExtra::grid.arrange(gg_cb, gg_ci, ncol = 1)