Continue predictor
Om het toxicologisch effect van een stof te kwantificeren worden dierproeven uitgevoerd waarbij dieren bij verschillende concentraties van de stof worden blootgesteld voor een bepaalde tijd. In dit voorbeeld wordt koolstofdisulfide (CS\(_2\)) bestudeerd in proeven met kevers. De centrale onderzoeksvraag is of de concentratie van CS\(_2\) een effect heeft op de mortaliteit (i.e. kans op sterven) van de kevers?
Design In 32 onafhankelijk experimenten wordt een kever blootgesteld
aan één van 8 concentraties (mg/l) van CS\(_2\) voor een gegeven
periode. De uitkomst van het experiment is: de kever sterft (\(y=1\)) of
de kever overleeft (\(y=0\)). Een R object met de data is opgeslagen in
de file kevers.rda in de dataset folder.
beetles <- read_csv("https://raw.githubusercontent.com/statomics/sbc20/master/data/beetles.csv")
head(beetles)
## # A tibble: 6 x 2
## dose status
## <dbl> <dbl>
## 1 169. 1
## 2 169. 0
## 3 169. 0
## 4 169. 0
## 5 172. 1
## 6 172. 0
table(beetles$dose, beetles$status)
##
## 0 1
## 169.07 3 1
## 172.42 3 1
## 175.52 3 1
## 178.42 2 2
## 181.13 1 3
## 183.69 0 4
## 186.1 0 4
## 188.39 0 4
We bouwen nu een logistisch regressiemodel waarbij we de log odds modelleren in functie van de dosis \(x_i\):
\[\log \frac{\pi_i}{1-\pi_i}=\beta_0+\beta_1 \times x_i.\]beetleModel<-glm(
status~dose,
data = beetles,
family = binomial)
summary(beetleModel)
##
## Call:
## glm(formula = status ~ dose, family = binomial, data = beetles)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7943 -0.7136 0.2825 0.5177 2.1670
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -53.1928 18.0046 -2.954 0.00313 **
## dose 0.3013 0.1014 2.972 0.00296 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 42.340 on 31 degrees of freedom
## Residual deviance: 26.796 on 30 degrees of freedom
## AIC: 30.796
##
## Number of Fisher Scoring iterations: 5
Het intercept heeft als betekenis de log odds op mortaliteit wanneer er geen \(\text{CS}_2\) gas wordt toegediend. Het duidt op een heel erg lage odds op sterfte (\(\pi/(1-\pi)=\exp(-53.2)\)) en dus op een kans die nagenoeg nul is. We verwachten inderdaad dat de kevers niet zullen sterven als we geen gas gebruiken in het experiment. Merk op dat de interpretatie van het intercept echter wel op een heel sterke extrapolatie berust gezien de minimum dosis in de dataset 169.07 mg/l bedroeg.
De geschatte odds ratio voor het effect van dosis op de mortaliteitskans is \(\exp(0.3013)=1.35\). Dus bij een toename van de dosis CS\(_2\) met 1 mg/l, is de odds ratio voor de mortaliteit \(1.35\). We besluiten dat dit effect heel significant is (\(p=\) 0.003). Een toename in de CS\(_2\) dosis doet de kans op sterven toenemen.
Figuur 59 toont de geschatte probabiliteit
in functie van de dosis die werd gefit a.d.h.v. het regressiemodel. De
figuur werd bekomen met onderstaande R code. We kunnen predicties voor
nieuwe dosissen berekenen d.m.v. de predict functie. Het argument
type="response" geeft aan dat we de predicties op de
probabiliteitsschaal wensen. Wanneer we dat niet specifiëren genereert
het model predictie op de log odds schaal.
beetlesTab<-table(beetles) %>% data.frame
data.frame(
grid = seq(
min(beetles$dose),
max(beetles$dose),.1)
) %>%
mutate(piHat = predict(beetleModel,
newdata = data.frame(dose=grid),
type = "response")
) %>%
ggplot(aes(grid, piHat))+
geom_line() +
xlab("dose") +
ylab("probability (dead)") +
geom_text(
aes(x = dose %>%
as.character %>%
as.double,
y = status %>%
as.character %>%
as.double,
label = Freq),
beetlesTab %>%
filter(status==0)
) +
geom_text(
aes(x = dose %>%
as.character %>%
as.double,
y = status %>%
as.character %>%
as.double,
label=Freq),
beetlesTab %>%
filter(status==1)
)
Figuur 59: Fit van het logistische regressiemodel voor de kevers data. Het aantal kevers die dood/levend was per dosis is weergegeven met een cijfer.