101_camilla.qmd

# Técnica não paramétrica e paramétrica para experimento com um fator

```{r, echo = FALSE, results = "asis"}
box::use(./config/setup)
setup$chapter_authors(c("Camilla Castellar",
                        "Walmes Marques Zeviani",
                        "Louise Larissa May De Mio"))
```

## Contextualização e descrição do experimento base para esta análise

As podridões em frutos resultam em redução de produtividade, uma vez que
frutos doentes não podem ser comercializados. Para entender melhor o
ciclo da doença, este estudo se concentrou em investigar diferentes
estágios de crescimento de frutos, visando orientar a aplicação de
produtos durante o período de maior suscetibilidade do fruto à infecçao do patógeno. O objetivo
principal foi estimar e comparar o período de incubação de uma podridão
que afeta macieiras em cinco estágios de crescimento dos frutos.

O experimento foi conduzido em condições controladas, usando frutos
colhidos no campo em cinco tamanhos diferentes. Após desinfestação, cada
fruto foi ferido na região equatorial e inoculado com
uma suspensão de conídios do patógeno sob um ferimento feito com agulha histológica. Os frutos foram então
distribuídos aleatoriamente em uma sala de incubação com condições
controladas e uniformes, seguindo um delineamento inteiramente
casualizado com cinco tratamentos (os tamanhos) e oito repetições (frutos)
(consulte a Figura @fig-esquema).

```{r}
#| label: fig-esquema
#| echo: false
#| fig-cap: "Ilustração do experimento onde frutos de 5 tamanhos são inoculados com o patógeno sob o ferimento feito na região equatorial."
knitr::include_graphics("101_camilla/tamanhos.png")
```

Posteriormente à inoculação, cada fruto foi avaliado diariamente ao
longo de 39 dias, período em que os frutos poderiam desenvolver sintomas
da doença. O dia em que os sintomas apareciam era registrado e
atribuído um status de 1 aos frutos sintomáticos na tabela de registro.
Para os frutos que não apresentavam sintomas, era registrado o status 0
e atribuído dia 39 (último dia de avaliação), indicando censura à
direita. Portanto, a tabela de dados deste estudo (Tabela
@tbl-conjuntodedados) continha apenas a data de ocorrência dos sintomas
para frutos sintomáticos ou a última data de avaliação para frutos sem
sintomas. Quando se faz necessário trabalhar com tabelas contendo o
histórico diário dos frutos, é preciso filtrar os dados para manter
apenas a data em que ocorreu a transição da condição sem sintoma para a
com sintoma.

```{r}
#| label: tbl-conjuntodedados
#| echo: false
#| tbl-cap: "Conjunto de dados para análise de sobrevivência exibindo os 3 primeiros registros para cada tamanho de fruto. A variável `day` é o dia de ocorrência do sintoma e a variável `status` é o indicador do tipo de censura: 0 é censura à direita e 1 é sem censura. Como a duração do experimento foi de 39 dias, a censura à direita foi na mesma data para todos os frutos."
tb <- read.csv2("101_camilla/sobrev.csv")

tb |>
    dplyr::group_by(tamanho) |>
    dplyr::slice_head(n = 3) |>
    knitr::kable(digits = c(1, 1, 1, 1),
                 align = c("cccc"),
                 row.names = FALSE)
```

<!-- WALMES -->

## Análise exploratória dos dados

Os dados deste experimento estão armazenados em um arquivo CSV contendo
3 colunas: `tamanho`, `day` e `status` que indicam o tamanho do fruto, o
dia de ocorrência do sintoma e o indicador de tipo de censura.

```{r}
#| message: false

# Carrega pacotes.
library(tidyverse)
library(multcomp)
library(emmeans)
library(survival)
library(survminer)
```

```{r}
#| message: false

# Importação dos dados do arquivo CSV.
dados  <-  readr::read_csv2("./101_camilla/sobrev.csv",
                            col_types = cols())
attr(dados, "spec") <- NULL
str(dados)

# Tabela com medidas resumo.
dados |>
    group_by(tamanho) |>
    rstatix::get_summary_stats(day, type = "common")
```

Conforme indicam as medidas resumo, são 8 unidades experimentais por
tratamento. No maior tamanho (`tam5`), o maior tempo para ocorrência dos
sintomas foi de 5 dias. Nos demais tamanho sempre teve ao menos um fruto
que levou 39 dias. Os frutos de maior tamanho apresentaram sintomas, no geral, antes que os frutos menores.

```{r}
# Ordena pelo dia do desfecho para fazer o gráfico ECDF.
dados0 <- dados |>
    group_by(tamanho) |>
    arrange(day) |>
    mutate(ord = seq_along(day)/n()) |>
    ungroup()
```

Para os tratamentos com tamanhos 1, 2 e 3 não se tem período de
incubação amostral pois, menos de 50% das repetições apresentaram
sintomas (Figura @fig-aed-frutos) no período de avaliação. No caso de poucos frutos
apresentarem sintomas, recomenda-se aumentar o número de repetições do
experimento. No entanto, isso não necessariamente impede de estimar a
curva de sobrevivência e estimar o período de incubação por meio de
modelos. A definição do período de incubação é o tempo necessário para
mais que metade dos indivíduos expostos apresentem o evento, ou seja, sintomas. Dessa forma, o
período de incubação é a mediana da variável aleatória tempo para a
ocorrência do evento. Com o ajuste do modelo aos dados, teremos a
estimativa para o período de incubação.

```{r}
#| label: fig-aed-frutos
#| fig-cap: "Análise exploratória do experimento. Cada linha horizontal é um fruto do experimento, ordenados conforme o tempo em que o sintoma foi observado. Os frutos estão agrupados pelo tratamento ao qual pertencem. Símbolos abertos indicam censura à direita e símbolos fechados indicam não censura."

# Gráfico com a linha do day de cada fruto que mostra quando o evento
# aconteceu e que tipo de desfecho.
ggplot(data = dados0,
       mapping = aes(x = day, y = ord, shape = factor(status))) +
    facet_grid(facets = tamanho ~ .) +
    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 = "Frequência acumulada") +
    scale_x_continuous(breaks = seq(0, max(dados$day), by = 10)) +
    theme(legend.position = c(0.95, 0.05),
          legend.justification = c(1, 0))
```

## Análise não paramétrica

Pela abordagem não paramétrica de análise de sobreviência, o período de
incubação de frutos como função dos tamanhos será obtido e comparado por
meio das curvas de Kaplan-Meier. O teste de log-rank será aplicado para
testar a hipótese de igualdade dos períodos de incubação.

Para começar, nos fragmentos a seguir, define-se o objetos com
diferentes representações de censura. No objeto `s1` definiu-se censura
à direita, o que não é exatamente o que ocorre no experimento. No objeto
`s2` a censura intervalar foi definida, que é o ocorre no experimento,
já que as avaliações são feitas 1 vez ao dia. É comum encontrar análises
usando censura à direita no lugar de intervalar quando a duração dos
intervalos é curta e de mesmo tamanho.

```{r}
# Resposta com avaliações sem censura e com censura à direita no final.
s1 <- with(dados,
           Surv(time = day,
                event = status))

# Resposta com avaliações com censura intervalar de 1 dia e com censura
# à direita no final.
s2 <- with(dados,
           Surv(time = day,
                time2 = day + 1,
                event = ifelse(status == 1, 3, 0),
                type = "interval"))

tibble("censura à direita (s1)" = s1, "censura intervalar (s2)" = s2)
```

O modelo pelo método de Kaplan-Meier foi ajustado aos dados usando a
resposta definida com a represetação de censura intervalar. Na saída
tem-se que a mediana não é estimada para 3 dos 5 tratamentos. O gráfico
com as curvas de Kaplan-Meier indicam que a sobrevida dos frutos maiores
é inferior aos demais.

```{r}
# Curvas sobreviência de Kaplan-Meier.
fit <- survfit(s2 ~ tamanho, data = dados)
fit

# Gráfico padrão com as curvas de Kaplan-Meier.
# plot(fit)

# Gráfico com melhor visual, legenda, etc.
ggsurvplot(fit,
           data = dados,
           # conf.int = TRUE,
           surv.median.line = "h",
           legend = "right", 
           ggtheme = theme_gray())

# Aplica o teste de log-rank. Só aceita censura à direita.
logrank <- survdiff(s1 ~ tamanho, data = dados)
logrank
```
```{r}
#| include: false
chi <- sprintf("%0.1f", logrank$chisq)
```

O teste do log-rank para testar a hipótese nula de igualdade entre as
curvas apresentou uma estatística de teste igual a `r chi`, bastante
desfavorável ($p < 0.001$) à hipótese nula, ou seja, em termos práticos
existe efeito do tamanho do fruto no comportamento da variável período
de incubação. Importante mencionar que o teste do log-rank não é feito
com dados de censura intervalar e sim os a censura à direita.

Apesar de ser possível testar o efeito de tratamentos, essa abordagem
não necessariamente irá estimar a mediana para todos os níveis, pois
ocorreu alta proporção de censura em alguns tamanhos de fruto neste
experimento. Como pode ser visto na saída da `survfit()`, apenas os
tamanhos 4 e 5 tiveram estimativa de mediana calculada.

A abordagem não paramétrica, conforme visto aqui, é mais recomendada
quando o número de unidades experimentais por tratamento for
suficientemente grande de forma que combinada com a prorporção de
censura, ainda permita estimar o tempo mediano. A censura pode ser
mitigada aumentando o tempo de avaliação do experimento para que se
observe o evento em mais unidades experimentais. Já que há dificuldade
de determinar a taxa de censura antes do experimento, as opções são:
usar um tamanho de amostra grande ou aumentar o período de monitoramento
do desfecho, sintoma no caso. O primeiro caso é uma decisão antes do experimento começar.
O último, pode ser feito a medida que o experimento é conduzido.

## Análise paramétrica

Na abordagem paramétrica vamos assumir uma distribuição de probabilidade
para o tempo de ocorrência dos sintomas. Existem várias distribuições
que podem ser consideradas, como Lognormal, Gama, Exponencial, Weibull,
entre outras. A distribuição Weibull é muito utilizada em análise de
sobrevivência e por isso será utilizada para exemplificar o uso da
abordagem paramétrica.

```{r}
#| eval: false

# Documentação disponível no R sobre os modelos paramétricos.
help(survreg.distributions, help_type = "html")
help(dweibull, help_type = "html")

survreg.distributions$weibull
```

A distribuição Weibull é uma distribuição de 2 parâmetros (existem
versões com mais parâmetros). A função de densidade da distribuição
usada na `dweibull()` e outras funções relacionadas é
$$
f(y) = \frac{a}{b} \left( \frac{y}{b} \right)^{a-1}
  \exp \left\{ -\left( \frac{y}{b} \right)^{a}\right\},
  \quad y > 0, a > 0, b > 0,
$$
em que $a$ é o parâmetro de forma (shape) e $b$ é o parâmetro de escala
(scale). Nessa parametrização, tem-se que:

  * Média: $\text{E}(Y) = b \Gamma (1 + 1/a)$.
  * Variância: $\text{Var}(Y) = b^2 (\Gamma(1 + 2/a) - (\Gamma(1 + 1/a))^2)$.
  * Função de distribuição: $F(y) = 1 - \exp\{-(x/b)^a\}$.

A parametrização usada pela `survreg()` é diferente mas matematicamente
relacionada a anterior por meio das expressões
$$
\begin{aligned}
  a &= 1/\lambda \quad &\Rightarrow \lambda &= 1/a \\
  b &= \exp\{\mu\} \quad &\Rightarrow \mu &= \log(b).
\end{aligned}
$$

Dessa maneira, o modelo estatístico para o tempo até aparecimento de
sintoma ($Y$) pode ser escrito como
$$
\begin{aligned}
  y_{ij} &= \text{Weibull}(a = 1/\lambda, b = \exp\{\mu_i\}) \\
  \mu_i &= \mu + \alpha_i \\
  \lambda &\propto 1,
\end{aligned}
$$
em que $\mu_i$ acomoda o efeito do nível $i$ de tamanho no parâmetro de
escala, $\lambda$ está vinculado ao parâmetro de forma. A não rejeição
da hipótese nula $H_0: \alpha_i = 0$ para todo $i$ corresponde não haver
efeito do tamanho dos frutos.

Nesta especificação de modelo, os parâmetros são estimados pelo método
da máxima verossimilhança. A função `survreg()` é usada para estimar o
modelo com distribuição Weibull para resposta com censura intervalar,
realizar o teste de hipótese e demais desdobramentos usuais de análise
de sobreviência.

Serão ajustados dois modelos para realizar o teste de hipótese para
efeito do tamanho de fruto. O modelo inicial não contém o efeito de
nunhuma fonte de variação e corresponde, portanto, o modelo nulo ou
modelo obtido sob a não rejeição da hipótese nula $H_0: \alpha_i = 0$
para todo $i$. O modelo seguinte acomoda o efeito de tamanho de frutos e
corresponde a especificação feita acima. O teste da razão de
verossimilhanças para esses modelos testa a hipótese nula por meio da
estatística diferença de deviances.

```{r}
#| message: false

# dados$day_orig <- dados$day
# dados$day <- log(dados$day_orig)

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

# Inclui o efeito de tamanho do fruto.
m1 <- update(m0, formula = . ~ tamanho)

# Valores de log-verossimilhança.
c(logLik(m0), logLik(m1))

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

Conforme indicado pelo teste da razão de verossimilhanças, existe
considerável evidência à favor da rejeição da hipótese nula, ou seja, a
favor do modelo que acomoda efeito de tamanho de fruto ($p < 0.001$). Em
termos práticos, isso quer dizer que existe efeito do tamanho de fruto
no tempo de incubação.

O fragmento de código abaixo tem propósito didático, pois mostra como
obter o valor de log-verossimilhança do modelo ajustado. Dessa forma,
fica registrado como estão relacionadas as parametrizações da Weibull na
`dweibull()` e `survreg()`, bem como a contribuição proveniente de dados
censurados e não censurados para a log-verossimilhança.

```{r}
# Avaliando a log-verossimilhança manualmente para o modelo `m0`.
a <- 1/m0$scale    # Parâmetro de forma da survreg() para a *weibull().
b <- exp(coef(m0)) # Idem pro parâmetro de escala.

# Log-verossimilhança e suas partes para censura intervalar.
sum(log(pweibull(dados$day,    # Parte dos dados com censura intervalar.
                 shape = a, scale = b,
                 lower.tail = FALSE) -
        pweibull(dados$day + 1,
                 shape = a, scale = b,
                 lower.tail = FALSE)) * dados$status) +
    sum(pweibull(dados$day,    # Parte dos dados censurados à direita.
                 shape = a, scale = b,
                 lower.tail = FALSE, log = TRUE) * (1 - dados$status))
```

Para estudar com mais detalhes o efeito do tamanho do fruto, serão
determinadas as estimativas dos parâmetros para obter a curva de
probabilidade de cada nível de tamanho. Com as estimativas individuais,
comparações múltiplas de parâmetros serão realizadas para testar
hipóteses de igualdade da distribuição do período de incubação entre
tamanhos de frutos.

As estimativas por tratamento para o parâmetro de escala são facilmente
obtidas com a função `emmean::emmeans()`. Tais estimativas não são
imeditamente interpretáveis, de modo que é mais interessante obter as
estimativas de média ou mediana a partir dos mesmos. Todavia, se houver
igualdade quanto ao parâmetro de escala, haverá para os demais, pois são
funções matemáticas uns dos outros. No código abaixo obtem-se a
estimativa de média.

```{r}
# Médias ajustadas (é na interpretação do parâmetro da Weibull).
emm <- emmeans(m1, specs = ~tamanho)
emm

# Estimativas conforme parametrização da {p, d, q, r}weibull.
tb_emm <- as.data.frame(emm)
a <- 1/m1$scale          # Shape.
b <- exp(tb_emm$emmean)  # Scale.

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

# Obtenção da mediana via função predict().
tb_emm$median_day <-
    predict(m1,
            type = "quantile",
            p = 0.5,
            newdata = tb_emm[, "tamanho", drop = FALSE])
tb_emm
```

```{r}
#| include: false

p1 <- sprintf("%0.2f", tb_emm$mean_day[5])
p2 <- sprintf("%0.2f", tb_emm$median_day[5])
```

Para os frutos de maior tamanho, o período de incubação médio é de `r
p1` dias e o período de incubação mediano é `r p2` dias. O período de
inbubação mediano corresponde à definição de período de inbubação
amostral no sentido de ser aquele para o qual metade das unidades
experimentais apresentam sintoma. Detalhe: com o atual parâmetro de
forma, a distribuição Weibull é assimétrica à direita e assim a mediana
é menor que a média.

Para fazer uma contemplação visual, serão obtidas as curvas de
probabilidade acumuladas. Essas curvas descrevem a probabilidade de
apresentar o sintoma em função do tempo saindo de uma probabilidade 0 no
instante 0 até uma probabilidade limite de 1 quando o tempo vai para
infinito. Pela análise da curva é possível determinar qual o tempo
necessário para que, por exemplo, 70% dos frutos apresentem sintoma. A
mediana, no caso, é o tempo para que 50% dos frutos apresentem sintoma.

```{r}
# Obtém as curvas estimadas da função de densidade e de distribuição.
tb_curves <- tb_emm |>
    group_by(tamanho) |>
    do({
        day <- seq(0, max(dados$day), length.out = 71)
        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 =  ~ tamanho) +
#     geom_line()

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

# Valores observados com sobreposição das curvas ajustadas.
ggplot(data = dados0,
       mapping = aes(x = day, y = ord, shape = factor(status))) +
    facet_grid(facets = tamanho ~ .) +
    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))
```

O gráfico exibe tanto os valores observados quanto as curvas estimadas
para cada tratamento. As curvas acompanham o padrão visto nos dados, ou
seja, as curvas de menor mediana são para os tratamentos em que os
frutos apresentaram sintomas mais cedo.

Com o intuito de melhor detalhar e efeito do tamanho do fruto, serão
testados os contrastes com o parâmetro de locação. Considerando todos os
contrastes par a par (contrastes de Tukey), são 10 hipóteses sob
investigação. A correção para a multiciplidade já é levada em conta na
implementação disponível na função `cld()`, o que torna muito fácil a
realização de procedimentos de comparação múltipla inclusive para
modelos que não sejam gaussianos.

O fragmento de código abaixo compara, portanto, os tamanhos de fruto em
relação ao parâmetro de escala. Os resultados são exibidos de forma
gráfica para mais rápida compreensão.

```{r}
# Comparando os tratamentos.
# Tabela com o resultado das comparações múltiplas.
tb_contr <- 
    emmeans(m1, specs = ~tamanho) |>
    multcomp::cld(Letters = letters) |>
    as.data.frame()
# tb_contr

# Puxa as demais estimativas para a mesma tabela.
tb_contr <- tb_contr |>
    rename("let" = ".group") |>
    inner_join(tb_emm[, c("tamanho", "mean_day", "median_day")]) |>
    mutate(let = trimws(let))
tb_contr |>
    arrange(tamanho) |>
    mutate_if(is.numeric, round, digits = 2)
```

```{r}
#| fig-height: 3

ggplot(data = tb_contr,
       mapping = aes(y = tamanho, x = emmean)) +
    geom_point() +
    geom_errorbarh(mapping = aes(xmin = lower.CL,
                                 xmax = upper.CL),
                   height = 0.1) +
    geom_text(mapping = aes(label = sprintf("%0.2f %s",
                                            emmean,
                                            let)),
              vjust = -0.5,
              hjust = 0.5) +
    labs(y = "Tamanho do fruto",
         x = "Estimativa do parâmetro de locação")
```

O gráfico mostra que não existe diferença entre os tamanhos 1, 2 e 3.
Para estes, a estimativa da mediana foi superior a 45 dias. O tamanho 5
foi o que apresentou menor tempo para os sintomas sendo diferente de
todos exceto do tamanho 4.

## Conclusão

  1. A abordagem não paramétrica:
     * A abordagem não paramétrica é interessante por não fazer muitas
       suposições sobre o fenômedo de interesse.
     * Ela permite facilmente testar o efeito global de uma fonte de
       variação.
     * No entanto, isso pode ser restritivo quando o experimento tiver
       mais fatores experimentais ou estruturais (e.g. blocos).
     * A abordagem não paramétrica não necessariamente irá estimar a
       mediana que é o período de incubação.
  2. A abordagem paramétrica:
     * Assume uma distribuição de probabilidade para o tempo entre
       eventos.
     * É capaz de acomodar os vários tipos de censura.
     * Permite acomodar o efeito de vários termos experimentais ou
       estruturais.
     * Permite estimar funções dos parâmetros, como média e mediana.
     * No entanto, se o modelo não for apropriado para o caso em
       questão, as inferências estatão comprometidas.

Nas análises aqui apresentadas, tem-se algumas considerações de itens
que podem ser avaliados mais a fundo.

O intuito da análise paramétrica foi demonstrar os passos da abordagem e
usou-se a distribuição Weibull. Todavia, outras distribuições poderiam
ser mais apropriados para estes dados. Caba ao pesquisador especificar e
ajustar diferentes distribuições e comparar a bondade de ajuste. Em
amostras pequenas, pode-se não se chegar a conclusões muito claras sobre
a melhor distribuição.

Uma forma de diagnósticar problemas em relação a distribuição assumida o
quanto a especificação do modelo é por meio da análise dos resíduos. A
análise dos resíduos e outras ferramentas de diagnóstico podem orientar
sobre mudar a distribuição ou aplicar uma transformação do tempo, como
logaritmica ou potência. Está fora do escopo entrar dos detalhes da
análise de diagnóstico para essa classe de modelo.