Hackeando o R

Hackeando o R: Dados censurados

By 3 de novembro de 2021 julho 28th, 2022 No Comments

No Hackeando o R de hoje, vamos apresentar métodos que são complementares ao post de semana passada. Enquanto a análise de sobrevivência era focada em observações que demoravam tempo demais para ocorrer (censuradas "à direita"), os métodos de hoje se interessam em lidar com o problema de dados que são censurados à esquerda. Isso pode ocorrer num contexto de testes de laboratório, como a mensuração da quantidade de certos elementos químicos na água, sendo que alguns dados podem não ser corretamente especificados por estarem em quantidade muito baixa, abaixo do limite de detecção das ferramentas do laboratório.

Os métodos apresentados no post sobre análise de sobrevivência conseguiriam (em tese) lidar com esse problema, porém na prática existem algumas diferenças que devem ser levadas em conta. Vamos começar nosso código carregando o pacote EnvStats, que contém funções construídas para lidar com dados não-detectados.

library(EnvStats)

Digamos que você possui observações de uma população que acredita que tem distribuição normal, porém os dados são censurados à esquerda. No código abaixo, geramos uma amostra nesse formato, criando 50 observações de uma normal com média 10 e desvio-padrão 3, e estabelecendo um valor mínimo de 8 para qualquer observação. Gravamos a média e desvio-padrão amostrais com e sem censura, para compararmos as estimativas posteriores.

# Normal(10, 3^2)
x = rnorm(5000, 10, 3)

mean_full = mean(x)
sd_full = sd(x)

# censura
censored = 8
x[censored] = 8

mean_cens_naive = mean(x)
sd_cens_naive = sd(x)

Para estimarmos a distribuição original, utilizaremos a função enormCensored, apropriada para a amostra dado que supomos que ela provém de uma dist. normal. Utilizaremos três métodos diferentes, cujos detalhes vão além do conteúdo desse post, e compararemos com o resultado de se estimar sem levar em conta a censura à esquerda. A função precisa apenas dos dados e de um vetor que indica quais observações foram censuradas.

estimates_mle = enormCensored(x, censored,
method = "mle")

estimates_bcmle = enormCensored(x, censored,
method = "bcmle")

estimates_ros = enormCensored(x, censored,
method = "ROS")

mean_full; sd_full
[1] 9.951524
[1] 2.975333
mean_cens_naive; sd_cens_naive
[1] 10.40674
[1] 2.337386
estimates_mle$parameters
 mean sd 
9.970963 2.942115 
estimates_bcmle$parameters
 mean sd 
9.970817 2.942756 
estimates_ros$parameters
 mean sd 
9.951850 2.968289 

Como podemos ver, as estimativas de todos os modelos ficam muito melhores do que a estimativa 'ingênua'. Abaixo, visualizamos o histograma dos dados, em conjunto com a distribuição normal estimada pelo modelo.


library(ggplot2)

ggplot(data.frame(x), aes(x)) +
geom_histogram(colour = "black", aes(y = ..density..)) +
geom_function(fun = dnorm, colour = "red", size = 2,
args = list(mean = estimates_mle$parameters[1], 
sd = estimates_mle$parameters[2]))

Quanto ao caso de estarmos interessados em uma regressão com dados censurados, podemos utilizar o modelo Tobit, apropriado para isso. Vamos utilizar novamente o pacote survreg, para analisar o dataset EPA.92c.zinc.df do pacote EnvStats. O dataset contém mensurações de zinco em 5 poços, sendo composto de 8 observações para cada um, que foram feitas (por hipótese) em momentos diferentes.


library(survival)

m = survreg(Surv(Zinc, !Censored, type = "left") ~ Sample + Well,
data = EPA.92c.zinc.df, dist = "gaussian")
summary(m)

Call:
survreg(formula = Surv(Zinc, !Censored, type = "left") ~ Sample +
Well, data = EPA.92c.zinc.df, dist = "gaussian")
Value Std. Error z p
(Intercept) 2.435 3.422 0.71 0.477
Sample2 4.854 3.508 1.38 0.166
Sample3 4.012 3.554 1.13 0.259
Sample4 5.872 3.495 1.68 0.093
Sample5 3.036 3.571 0.85 0.395
Sample6 4.764 3.513 1.36 0.175
Sample7 5.655 3.501 1.62 0.106
Sample8 2.345 3.598 0.65 0.514
Well2 1.474 2.578 0.57 0.568
Well3 1.388 2.578 0.54 0.590
Well4 2.632 2.542 1.04 0.301
Well5 1.135 2.578 0.44 0.660
Log(scale) 1.497 0.181 8.29 <2e-16

Scale= 4.47

Gaussian distribution
Loglik(model)= -72 Loglik(intercept only)= -74.7
Chisq= 5.58 on 11 degrees of freedom, p= 0.9
Number of Newton-Raphson Iterations: 3
n= 40

________________________
(*) Para entender mais sobre análises estatísticas, confira nosso Curso de Estatística usando R e Python.

Receba diretamente em seu e-mail gratuitamente nossas promoções especiais
e conteúdos exclusivos sobre Análise de Dados!

Assinar Gratuitamente