A toxicologia estuda os efeitos prejudiciais de substâncias químicas em organismos. Diversos testes de toxicidade podem ser usados avaliar o impacto dessas substâncias nos organismos [1]. Quando se analisa a mortalidade em testes de toxicidade, é comum o cálculo da concentração letal para 50% da população testada, conhecida como CL50. Em alguns casos, também é importante determinar outras concentrações, a exemplo da CL10 ou CL90. Uma ferramenta estatísitca comumente usada neste tipo de análise é o modelo de regressão de Probit [2].
Neste texto, uma análise de Probit será realizada para calcular concentrações letais em linguagem R. Um conjunto de dados da literatura, disponível com livre acesso na internet, será usado como caso de estudo. Para esta aplicação, primeiramente serão apresentados conceitos básicos sobre modelos lineares generalizados e o modelo de Probit, além das distribuições de probabilidade importantes para a análise.
Os modelos lineares generalizados (GLM) são extensões flexíveis dos modelos de regressão linear, que permitem que a variável resposta possua uma distribuição de erro diferente da normal. Os GLM são formados por um componente aleatório, um componente sistemático e uma função de ligação:
O modelo de regressão de Probit é um caso especial de GLM, que usa como função de ligação a função de distribuição normal padrão \(\Phi(\cdot)\), conforme mostrado a seguir.
\[\mu_i=\Phi(\eta_i)=\Phi(\beta_0 + \beta_1X_{1i} + \dots + \beta_kX_{ki} + u)\]
No modelo Probit, o valor do preditor linear \(\eta_i\) é tido como um quantil \(z\) da distribuição normal padrão, isto é, \(\Phi(z)=P(Z\leq z)\), em que \(Z \sim Normal(0,1)\). A variável resposta é convertida em uma variável contínua que assume valores reais no intervalo \([0,1]\) por meio da função de ligação. Assim, o modelo Probit pode ser pensado como um modelo linear para a transformação da resposta ou como um modelo de regressão não linear para a resposta [4].
Seja \(Y_i\) uma variável binária, cujas realizações podem assumir somente dois valores possíveis (sucesso se \(Y_i=1\) ou fracasso se \(Y_i=0\)), pode-se usar a distribuição de Bernoulli para representá-la. Então, diz-se que \(Y_i\sim Be(p)\) com média \(\mu_i=p\), em que \(p\) é a probabilidade de sucesso. Na distribuição de Bernoulli, o experimento é realizado apenas uma vez. No caso em que o experimento seja realizado mais de uma vez, assumindo que as realizações são independentes entre si, pode-se usar a distribuição Binomial. Esta distribuição descreve o comportamento dos resultados de um número \(n\) de realizações de um experimento com probabilidade de sucesso \(p\). Neste caso, diz-se que \(Y_i\sim Bin(n,p)\) com média \(\mu_i=np\) [5].
Os resultados de testes de concentração letal são a relação de organismos mortos por diferentes concentrações de alguma substância tóxica a qual foram expostos. Pode-se representar esta relação por duas variáveis binárias, que indicam respectivamente a morte e a sobrevivência de cada organismo da amostra. Por exemplo, para um organismo, tem-se dead = 1 e alive = 0 caso haja a morte, e vice-versa. Em vez de visualizar o resultado individualmente, uma outra opção é considerar o resultado de todos os organismos, somando-se o total obtido em cada variável. Neste caso, o resultado será dado por duas variáveis que indicam respectivamente o total que sobreviveu ou não. Por exemplo, em uma amostra formada por 10 organismos, tem-se dead = 6 e alive = 4. Assim, as variáveis binárias são convertidas em variáveis de uma distribuição Binomial.
Neste texto, foi usado um conjunto de dados da Mendeley, disponíveis em [6]. São dados de testes de toxicidade de concentração letal (CL50) realizados com espécies de nematoides em resposta ao cloreto de cobre e cloreto de zinco. Os arquivos contém dados brutos de toxicidade, cálculos de CL50 e validação de concentração química. Em particular, neste texto, foram usados os dados de cepas N2 de C. elegans submetida a cloreto de zinco.
No R, foi usado uma variável para o número de mortes (dead) e outra para o número de sobreviventes (alive), além da concentração de zinco (conc), medida em mg/L. Em seguida, o conjunto de dados experiment foi criado como um data.frame.
# Resultados do experimento
dead <- c(0, 0, 0, 0, 3, 3, 2, 2, 5, 5, 4, 5, 3, 9, 6, 6, 8, 10, 6, 7, 9, 8, 9, 8, 8, 9, 9, 9)
alive <- c(10, 10, 10, 10, 7, 7, 8, 8, 5, 5, 6, 5, 7, 1, 4, 4, 2, 0, 4, 3, 1, 2, 1, 2, 2, 1, 1, 1)
conc <- c(0, 0, 0, 0, 34.5744, 34.5744, 34.5744, 34.5744, 68.5008, 68.5008, 68.5008, 68.5008, 136.3536, 136.3536, 136.3536, 136.3536, 204.2064, 204.2064, 204.2064, 204.2064, 272.0592, 272.0592, 272.0592, 272.0592, 339.912, 339.912, 339.912, 339.912)
# Cria o conjunto de dados
experiment <- data.frame(dead, alive, conc)## dead alive conc
## 0 10 0.0000
## 0 10 0.0000
## 0 10 0.0000
## 0 10 0.0000
## 3 7 34.5744
## 3 7 34.5744
## 2 8 34.5744
## 2 8 34.5744
## 5 5 68.5008
## 5 5 68.5008
## 4 6 68.5008
## 5 5 68.5008
## 3 7 136.3536
## 9 1 136.3536
## 6 4 136.3536
## 6 4 136.3536
## 8 2 204.2064
## 10 0 204.2064
## 6 4 204.2064
## 7 3 204.2064
## 9 1 272.0592
## 8 2 272.0592
## 9 1 272.0592
## 8 2 272.0592
## 8 2 339.9120
## 9 1 339.9120
## 9 1 339.9120
## 9 1 339.9120O modelo Probit foi construído usando a função glm do próprio R. Nesta função, o argumento family especifica a distribuição de probabilidades binomial com função de ligação (link function) probit.
# Constrói o modelo Probit
mod <- glm(
formula = cbind(dead, alive) ~ conc,
data = experiment,
family = binomial(link = "probit")
)Neste texto, foi usada a formulação de Wilkinson-Rogers [7] para construção do GLM, que usa a notação cbind(dead, alive) para definir as entradas do modelo. Em R, existem outras alternativas para a definição das entradas da função glm.
Construído o modelo, pode-se verificar uma breve descrição usando a função summary.
summary(mod)##
## Call:
## glm(formula = cbind(dead, alive) ~ conc, family = binomial(link = "probit"),
## data = experiment)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9373 -0.5917 -0.2195 0.4735 2.5891
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.9498449 0.1384382 -6.861 6.83e-12 ***
## conc 0.0074363 0.0008217 9.050 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 142.39 on 27 degrees of freedom
## Residual deviance: 40.30 on 26 degrees of freedom
## AIC: 100.09
##
## Number of Fisher Scoring iterations: 5Os coeficientes do modelo podem ser obtidos da seguinte maneira:
b0 <- mod$coefficients[1] # Intercept (b0)
b1 <- mod$coefficients[2] # Slope (b1)Para o cálculo da concentração letal, deve-se definir qual a porcentagem p que se deseja calcular. Por exemplo, para a CL50, tem-se p = 50. Em seguida, calcula-se o quantil da distribuição normal referente ao valor p usando a função qnorm. Por fim, o valor da concentração letal é calculado usando os coeficiente do modelo.
# CL50
p <- 50
# Calcula o quantil referente a CL50
est <- qnorm(p / 100)
# Calcula a CL50
CL50 <- (est - b0) / b1
# Valor da cL50 calculada
CL50## (Intercept)
## 127.7306Através do modelo construído, obteve-se um valor de CL50 igual a 127,7306 mg/L.
Os autores do conjunto de dados usado neste texto [6] optaram por um modelo log-logit para calcular a CL50 usando a função glm do pacote MASS. O resultado obtido por eles neste caso foi igual a 124,24 mg/L. Comparando as duas abordagens, os valores calculados foram próximos e a diferença observada pode ser devida ao modelo diferente usado pelos autores.
Este texto demostrou a aplicação de análise de Probit em linguagem R para o cálculo de concentrações letais em testes de toxicidade. Um conjunto de dados da literatura foi usado como caso de estudo e a análise de Probit foi usada com sucesso para estimar a concentração letal CL50 de zinco em uma cepa de uma espécie de nematoides. A diferença entre os resultados obtidos pelos autores do conjunto de dados pode ser devida ao modelo diferente usado por eles.