102-thiago.Rmd

# Análise de sobrevivência por abordagem paramétrica com experimento fatorial: caso de experiemntospara resistência prática com fungicidas.

```{r, echo = FALSE, results = "asis"}
chapter_authors(c("Thiago de Aguiar Carraro",
                  # "Paulo dos Santos Faria Lichtemberg",
                  "Walmes Marques Zeviani",
                  "Louise Larissa May De Mio"))
```

A análise de sobrevivência infere sobre a probabilidade de um evento
ocorrer no tempo. Neste capítulo será abordada a análise de experimentos
com mais de um fator. Dada a limitação em acomodar estruturas
experimentais visto na abordagem com Kaplam-Meier, a análise com modelos
paramétricos será utilizada.

Os modelos paramétricos podem estimar a distribuição do tempo de vida de
um patógeno e, ao mesmo tempo determinar parâmetros importantes para
compreender a interações de fatores [@Zhang2016]. Além do Weibull, visto
anteriormente, outros modelos de ajustes paramétricos podem ser
utilizados dependendo dos dados e/ou experimentos, como as regressões
exponenciais, gamma e log-logístico [@Zhang2016]. Portanto, este modelo
que será usado para análise dos dados de uma situação real (experimento
*ex vivo* com fungicidas) que será apresentado na sequência como caso 2
de análise de sobrevivência.

## Contextualização e descrição do experimento

O controle químico é uma das formas mais eficientes para o controle de
doenças de plantas, porém alguns fungicidas de sítio-específicos
apresentam riscos de médio a alto, como os IDMs (Inibidores da
DesMetilação) e IQes (Inibidores da Quinona externa), respectivamente,
para selecionar isolados do patógeno com menor sensibilidade. Esta
seleção pode levar a uma perda da eficiência de controle destes de
produtos comerciais com esses ingredientes ativos em campo. Para adotar
estratégias de controle eficientes, garantindo o manejo adequado ao
longo das safras, é necessário que seja realizado o monitoramento da
sensibilidade do fungo a fungicidas [@lichtemberg2016].

Para estudos de monitoramento da sensibilidade de fungo à fungicidas,
são utilizados ensaios *in vitro*, como: i) dose discriminatória (única
dose de fungicida capaz de separar os isolados em fenótipos sensiveís ou
resistentes), ii) CE~50~ (concentração efetiva que inibe 50% do
crescimento ou germinação do esporo do patógeno) e iii) CMI
- concentração mínima inibitória. Com esses dados, é possível
comparar o nível de sensibilidade dos isolados ao longo do tempo
(anos/safras) e verificar se está ocorrendo mudanças da sensbilidade do
patógeno aos fungicidas [@russell], para melhor adoção das estratégias
de controle das doenças.

Baseado nos ensaios *in vitro*, isolados são classificados
fenotipicamente em relação a sensibilidade aos fungicidas, e alguns
isolados representativos são selecionados para realização de ensaios
ensaios *ex vivo*, contrastando fenótipos em relaçao a sensibilidade a fungicidas.
Estes testes são recomendados para validar uma
possível resistência prática de fungicidas no controle de uma doença.
Resistência prática, portanto, está relacionada a uma possível falha do
controle da doença em condições de campo [@ghinirekimatih2000]. Estes
ensaios podem ser realizados com plantas (*in vivo*) ou partes da planta
destacadas (*ex vivo*), que serão submetidos a um tratamento com a
concentração recomendada na bula do produto para aplicação na área de
produção [@Schnabel2015].

Além das técnicas mais usuais nos trabalhos de monitoramento, a
observação do período de incubação (PI) em um ensaio de eficiência de
fungicidas pode ser uma medida epidemiológica importante para se
analisar o efeito de controle do produto ao longo do tempo, podendo
considerar este tempo o intervalo de uma aplicação à outra. Ou seja, o
PI vai fornecer informações do tempo efetivo de controle de um fungicida
a um patógeno, e isso é importante considerando isolados de regiões ou
safras diferentes. Uma vez detectada a mudança de sensibilidade em
populações de patógenos, novas recomendações e estratégias para o uso de
fungicidas devem ser adotadas.

A análise de sobrevivência proposta aqui, teve o o objetivo de avaliar
diferentes tratamentos (fungicidas) para o controle de duas espécies de
*Colletotrichum*, agentes causais da antracnose do caquizeiro no Brasil.
O experimento foi conduzido em esquema fatorial duplo $6 \times 2$
combinando por completo, sendo seis fungicidas (A, B, C, D, E e testemunha) com
duas espécies patogênicas (X e Y). Cada combinação experimental contou
com 25 repetições e o experimento foi realizado duas vezes.
Adequaçoes a depender dos fenótipos a serem testados, espécies ou tipos de fungicidas podem ser feitas para ajustar o delineamento e a análise.
Os frutos foram tratados com os produtos previamente, e após 24 horas
foram inoculados com as espécies fúngicas. Durante um total de 10 dias,
diariamente foi examinado se existia a presença ou não de lesão
(sintomas) nos frutos. Em caso positivo, era anotado a status um (1),
caso contrário era marcado como status zero (0). Sendo assim, o período
de incubação (PI) foi estimado quando metade ou mais das unidades
experimentais de uma combinação experimental (13 frutos, do total de 25 frutos) apresentaram o
sintoma da doença (Figura \@ref(fig:prancha-1)).

(ref:prancha-1) Esquema ilustrativo do ensaio de eficiência de fungicidas: os frutos foram previamente desinfestados com álcool (70%), hipoclorito (0,5%) e água destilada autoclavada, em seguida, os frutos foram tratados com fungicidas, e após 24 horas foram inoculados (com suspensão de 1x10^4^ conídios.mL^-1^). Na figura, as letras de A a E estão representando os tratamentos fungicidas (com 25 repetições) e a testemunha pela sigla 'test', enquanto as letras X e Y são os patógenos.

```{r prancha-1, echo = FALSE, fig.cap= '(ref:prancha-1)'}
knitr::include_graphics("./102-thiago/prancha_1.jpg")
```

## Análise de sobreviência

Nessa seção, será apresentada a abordagem paramétrica, assumindo a
distribuição de Weibull, para analisar o experimento com estrutura
fatorial.

Os dados estão armazenados em um arquivo de texto pleno (`*.txt`)
contendo 5 colunas. Após a leitura, categorizou os dados estabelecendo
os dias em que o evento ocorreu. Lembrando, sendo o evento determinado
pelo tempo onde 50% ou mais das repetições apresentaram sintomas da
doença para cada condição experimental. Em seguida, um gráfico
exploratório foi realizado para exibir esse tempo (dias) e o tipo de
desfecho de cada fruto para cada tratamento. A análise exploratória
sugere que tratamentos ineficientes no controle da doença apresentaram
tempo de incubação igual ou inferior a testemunha. Além disso, avalia a
resposta de duas diferentes espécies agentes causais em causar doença
quando inoculadas em cada fruto tratado.


```{r, message = FALSE}
# Pacotes.
rm(list = objects())

library(tidyverse)
library(survival)
library(emmeans)
```

O bloco de código abaixo faz a importação dos dados, conversão de tipo
de valor nas variáveis que são termos experimentais e mostra como estão
os registros para uma unidade experimental.

```{r, message = FALSE}
# Importação.
da <- read_tsv("102-thiago/frutosmaduros.txt")
attr(da, "spec") <- NULL
str(da)

# Tabela de frequência cruzada.
xtabs(~spp + fung, data = da)

# Converter termos experimentais em fatores.
da <- da %>%
    mutate(fung = factor(fung),
           spp = factor(spp),
           ue = interaction(fung, spp, rep))

# A tabela contém o registro dia após dia de cada unidade experimental.
# O evento ocorre quando o variável `status` troca de valor 0 para 1.
da %>%
    filter(ue == levels(ue)[1])
```

Como mostrado, a primeira UE, o fruto de rótulo `A.X.1`, teve o desfecho
observado no dia 2. Depois dessa data, os valores serão todos iguais a
1, bem como antes dela, todos iguais a 0. Essa forma de armezenar mostra
o status para cada dia de avaliação. Para análise de sobreviência só é
necessário o dia em que o status mudou de 0 para 1.

No bloco de código a seguir, os dados no formato longo com anotações
diárias da condição das UEs são transformados para a representação
necessária para a análise de sobrevivência.

Na prática, não é necessário anotar diariamente 0, 0, 0, etc, para datas
quando a unidade experimental (UE) ainda não apresentou o desfecho e 1,
1, 1, etc para datas após a UE apresentar o desfecho. Basta registrar as
estatísticas suficientes que são: identificador da UE ou condição
experimental e a data de ocorrência do evento desfecho. Essa data
corresponde ao dia em que observou-se aparecimento de sintoma no fruto,
neste caso. Aquelas UEs que não apresentaram sintomas até o final do
período experimental, são as UEs com resposta censurada, em particular,
censura à direita, visto que o evento irá certamente acontecer e ser
obervado se o período experimental for prolonagado.

Diante do que foi dito, o pesquisador pode optar por registrar os dados
no formato longo com anotações diárias por UE, que é o que temos, ou
anotar apenas as estatísticas necessárias, com um linha na tabela por
UE. O código abaixo mostra como passar do primeiro formato para o
segundo.

```{r, message = FALSE}
# Determina a primeira data de valor `status == 1` para cada unidade
# experimental.
da_1 <- da %>%
    group_by(ue) %>%
    filter(status == 1) %>%
    top_n(day, n = -1) %>%
    ungroup()
da_1

# Determina a última data de valor `status == 0` para cada unidade
# experimental restante, ou seja, as que não apresentaram `status == 1`.
da_0 <- anti_join(da, da_1, by = "ue") %>%
    group_by(ue) %>%
    top_n(day, n = 1) %>%
    ungroup()
da_0

# Junta as duas porções.
da_final <- bind_rows(da_1, da_0)
da_final

# Desenho experimental (5 + 1) * 2 com 25 frutos por condição
# experimental.
xtabs(~fung + spp, data = da_final)

# Número de frutos com desfecho observado.
xtabs(status ~ fung + spp, data = da_final)
```

Agora com a tabela contendo os dados necessários e no formato apropriado
para a análise de sobreviência, pode-se realizar a análise exploratória
para uma inspeção preliminar dos dados. É importante realizar essa etapa
para verificar eventuais inconsistências nos valores, erros de
preenchimento, perdas de unidiades experimentais, etc, bem como
antecipar alguns resultados que poderão ser confirmados pela análise.

```{r}
# Ordena pelo dia do evento dentro de cada unidade experimental.
da_final <- da_final %>%
    group_by(spp, fung) %>%
    arrange(day) %>%
    mutate(ord = seq_along(day)/n()) %>%
    ungroup()
```

(ref:ecdf-frutos) Distribuição acumulada empírica do tempo até ocorrência do desfecho nos frutos do experimento separados por condição experimental. Pontos indicam o tipo de desfecho (com censura ou não). As linhas vão desde o início do experimento até o momento de registro do desfecho. A linha horizontal está na posição 0.5 e permite determinar a mediada no tempo cuja interpretação é o período de incubação.

```{r ecdf-frutos, echo = FALSE, fig.cap= '(ref:ecdf-frutos)'}
# Gráfico com a linha do tempo de cada fruto que mostra quando o evento
# aconteceu e que tipo de desfecho.
ggplot(data = da_final,
       mapping = aes(x = day, y = ord, shape = factor(status))) +
    facet_grid(facets = spp ~ fung) +
    geom_point() +
    geom_segment(mapping = aes(x = 0, y = ord, xend = day, yend = ord)) +
    scale_shape_manual(breaks = c(0, 1),
                       values = c(1, 19),
                       labels = c("Censura", "Sintoma")) +
    labs(shape = "Desfecho",
         x = "Dias após inoculação",
         y = "Ordem do fruto") +
    scale_x_continuous(breaks = seq(0, max(da_final$day), by = 2)) +
    geom_hline(yintercept = 0.5, color = "orange")
```

Como são vários frutos por combinação experimental e o período de
observação amplo o suficiente, sempre obteve-se frutos apresentando o
desfecho. Houveram condições experimentais, no entanto, com mais da
metade dos frutos apresentando censura à direita.

Será feita análise de sobrevivência com esses dados, pois é motivada
pelo tipo de resposta de natureza "tempo até desfecho". Para acomodar a
estrutura de planejamento fatorial completo 2 x 6, será usada uma
abordagem paramétrica. Por apresentar boa flexibilidade, será usada a
distribuição Weibull para o tempo até o desfecho. Com isso consegue-se
acomodar a censura e a estrutura experimental.

```{r, eval = FALSE, include = FALSE, message = FALSE}
# Antes de ir para o ajuste da Weibull, vamos reconhecer a
# parametrização.

# Gera números aleatórios de uma distribuição Weibull.
y <- rweibull(n = 1000, shape = 4, scale = 2)
plot(ecdf(y))

# Ajusta o modelo com a `survreg()`.
m <- survreg(formula = Surv(time = y) ~ 1, dist = "weibull")
summary(m)

# Compreende os coeficientes estimados.
exp(coef(m)) # scale = exp(coef) -> 2.

1/m$scale    # shape = 1/m$scale -> 4.

# Verifica o ajuste com obsevado vs ajustado.
plot(ecdf(y))
curve(pweibull(x, shape = 1/m$scale, scale = exp(coef(m))),
      add = TRUE,
      col = 2)

# Obtenção do tempo médio.
a <- 1/m$scale
b <- unname(exp(coef(m))[1])
c("sample_mean" = mean(y),
  "estimated_mean" = b * gamma(1 + 1/a))
```

Nota-se que a análise de sobrevivência é a abordagem mais adequada para
análise destes dados, visto que de fato a variável resposta é o tempo
para o desfecho e tem-se censura. A modelagem para a variável aleotória
é o 'tempo para o aparecimento de sintomas (desfecho)'

Para testar a significância dos termos experimentais, uma sequência de
modelos aninhados será ajustada aos dados. O teste da razão de
verossimilhanças será aplicado para testar a contribuição de cada termo
experimental.

```{r, message=FALSE}
# Definição do tipo de censura.
s <- with(da_final,
          Surv(time = day,
               event = status))

# Modelo nulo que não tem efeito de nenhum fator.
m0 <- survreg(formula = s ~ 1,
              data = da_final,
              dist = "weibull")

# Inclui o efeito de espécie.
m1 <- survreg(formula = s ~ spp,
              data = da_final,
              dist = "weibull")

# Inclui o efeito de espécie e fungicidas (aditivo).
m2 <- survreg(formula = s ~ spp + fung,
              data = da_final,
              dist = "weibull")

# Inclui o efeito de espécie e fungicidas com interação.
m3 <- survreg(formula = s ~ spp * fung,
              data = da_final,
              dist = "weibull")

# Teste da razão de verossimilhanças entre modelos encaixados.
anova(m0, m1, m2, m3)
```

O modelo com interação é o mais apropriado para descrever o tempo para o
aparecimento de sintomas. O teste da razão de verossimilhanças entre
modelos aninhados indicou que o termo de interação é significativo.
Portanto, existe interação entre espécies e produtos ao nível de 5% pelo
teste de deviance entre modelos encaixados

```{r, message=FALSE}
# Modelo final.
mx <- m3

# Extração dos tempos médios e confecção das curvas de sobrevivência.
# Médias ajustadas (é na interpretação do parâmetro da Weibull).
emm <- emmeans(mx,
               specs = ~spp + fung)
emm

tb_emm <- as.data.frame(emm)

# Estimativas conforme parametrização da {p, d, q, r}weibull.
a <- 1/mx$scale          # shape.
b <- exp(tb_emm$emmean)  # scale.

# Obtenção do tempo médio (é uma função complexa dos dois parâmetros).
tb_emm$mean_day <- b * gamma(1 + 1/a)
tb_emm

# Curvas de probabilidade de desfecho.
tb_curves <- tb_emm %>%
  group_by(spp, fung) %>%
  do({
    day <- seq(0, 20, length.out = 51)
    dens <- dweibull(day, shape = a, scale = exp(.$emmean))
    accu <- pweibull(day, shape = a, scale = exp(.$emmean))
    data.frame(day, dens, accu)
  }) %>%
  ungroup()

# # As curvas de densidade ajustadas.
# ggplot(data = tb_curves,
#        mapping = aes(x = day, y = dens)) +
#   facet_grid(facets = spp ~ fung) +
#   geom_line()

# # As curvas de "1 - sobreviência".
# ggplot(data = tb_curves,
#        mapping = aes(x = day, y = accu)) +
#   facet_grid(facets = spp ~ fung) +
#   geom_line()

# Valores observados com sobreposição das curvas ajustadas.
ggplot(data = da_final,
       mapping = aes(x = day, y = ord, shape = factor(status))) +
  facet_grid(facets = spp ~ fung) +
  geom_point() +
  scale_shape_manual(breaks = c(0, 1),
                     values = c(1, 19),
                     labels = c("Censura", "Sintoma")) +
  labs(shape = "Desfecho",
       x = "Dias após inoculação",
       y = "Frequência acumulada") +
  geom_vline(data = tb_emm,
             mapping = aes(xintercept = mean_day),
             color = "orange") +
  geom_line(data = tb_curves,
            inherit.aes = FALSE,
            mapping = aes(x = day, y = accu))
```

Para finalizar a análise será realizado o estudo de interação dos
produtos em cada espécie inoculada por meio do desdobramento com testes
para o parâmetro de locação. Com esses valores podem fazer inferências
sobre o tratamento mais ou menos eficientes no controle da antracnose do
caquizeiro. Observa-se que existe uma diferença no controle da doença
entre as espécies.

(fig:errorbar-fungxesp) Estimativas para o parâmetro de locação da distribuição Weibull, acompanhadas do intervalo de confiança (95%), para cada uma das combinações experimentais.

```{r errorbar-fungxesp, message=FALSE, fig.cap = '(fig:errorbar-fungxesp)'}
# Estudar produtos em cada espécie.
emm_fung_in <- emmeans(mx, specs = ~fung | spp)
emm_fung_in

# contrast(emm_fung_in, method = "pairwise")
fung_in_comp <- multcomp::cld(emm_fung_in, Letters = letters)
fung_in_comp

# Tabela com o resultado das comparações múltiplas.
tb_contr <- as.data.frame(fung_in_comp) %>%
    mutate_at(".group", trimws)
tb_contr

ggplot(data = tb_contr,
       mapping = aes(x = fung, y = emmean)) +
  facet_wrap(facets = ~spp) +
  geom_point() +
  geom_errorbar(mapping = aes(ymin = lower.CL,
                              ymax = upper.CL),
                width = 0.05) +
  geom_text(mapping = aes(label = sprintf("%0.2f %s",
                                          emmean,
                                          .group)),
            hjust = 0.5,
            vjust = -1,
            angle = 90,
            nudge_x = 0.075) +
  labs(x = "fungicidas",
       y = "Estimativa do parâmetro de locação")
```

Para a espécie X, houve diferença da testemunha apenas com os fungicidas
B, C e E. Para a espécie Y, a houve diferença da testemunha apenas para
com o fungicida C.

Para estimar os temos médios e medianos, pode-se usar as instruções
apresentadas no capítulo anterior.

## Conclusão

A análise com a abordagem paramétrica é interessante para este tipo de estudo, além de ser
apropriada para o tipo de variável resposta em questão, é capaz de
acomodar a estrutura experimental descrita e outras combinações são possíveis.