Download Notes on Binomial Data - Examples with Resolution | MATH 6080 and more Study notes Mathematics in PDF only on Docsity!
Binomial Data
1 Binomial Data
1.1 Working example - Tobacco Budworms
Example 1. The mortality of the moth tobacco budworm “Heliothis virescens” was analyzed in an experiment, where groups of 20 individuals were exposed to the pyrethroid trans-cypermethrin and the number of deaths recorded after 72 hours. Six different doses were analyzed for each group of 20 male and female moths. The data are summarized as follows:
Mortality of moths exposed to cypermethrin.
Dose [μg] of cypermethrin Dead Males Dead Females 1 1 0 2 4 2 4 9 6 8 13 10 16 18 12 32 20 16
Note:
- Doses are doubled.
- Data suggests that the probability of dying increases with dose.
- Data suggests that males are more susceptible than females.
1.2 Binomial distribution
Example 2. A total of n = 5 female moths were exposed to a pesticide and it was recorded that y = 2 of these died.
Note:
- The data consists of only one observation since the aggregated “fate” of the moths (number of dead) is reported and not the fate of each single moth.
- To model this data we will make the following assumptions:
- The event of death is observed for a prefixed number of moths.
- Death of moths occur independently.
- Each moth has the same probability of dying denoted with π.
- Under these assumptions the number y of dead moths can be described using a bino- mial distribution, written y ∼ Bin(n, π).
Calculating probabilities using the binomial distribution.
Question: What is the probability of observing 2 dead moths (written P (y = 2))?
Let zi, for i = 1, 2 , 3 , 4 , 5, represent each of the five moths, respectively. The two dead moths out of the 5 could have occurred in the following 10 ways.
z 1 1 1 1 1 0 0 0 0 0 0 z 2 1 0 0 0 1 1 1 0 0 0 z 3 0 1 0 0 1 0 0 1 1 0 z 4 0 0 1 0 0 1 0 1 0 1 z 5 0 0 0 1 0 0 1 0 1 1 y =
∑ zi 2 2 2 2 2 2 2 2 2 2
Letting n = 5, y = 2, these 10 possibilities can be computed using combinations as follows: ( n y
)
n! y!(n − y)!
( 5 2
)
Note:
- The probability of death for moths 1 and 3 and survival of the others is given by
P (z 1 = 1, z 2 = 0, z 3 = 1, z 4 = 0, z 5 = 0) = π(1 − π)π(1 − π)(1 − π) = π^2 (1 − π)^3.
- The probability of the death of two arbitrary moths is
P (y = 2) =
( 5 2
) π^2 (1 − π)^3.
- An estimate for π is the observed proportion
πˆ =
y n
- Putting these together, an estimate for the probability of observing 2 dead moths out of 5 exposed moths is given by ( 5 2
) (.4)^2 (1 − .4)^3 =. 35.
- One can easily solve for the probability π from the η = logit(π) via the logistic function (based on the cdf of the logistic random variable) as follows:
π =
exp η 1 + exp η
(See GLM Overview notes for details on η = x′β. )
Example - In the experiment with n = 5 moths and y = 2 deaths, we have three different descriptions related to the mortality probability π given by:
Probability, odds, and logit of π.
ˆπ oddŝ ηˆ = logit ˆπ .4 .67 = (^5) −^22 -.
1.4 Comparison of Two Treatments
Example - Recall in the tobacco budworm experiment that moths were exposed to different doses of a poison by gender. If we consider just the male budworms which received doses 1 and 2 at which y 1 = 1 and y 2 = 4, respectively, out of n = 20 died, there are several options to compare the lethality of the doses.
- We assume that y 1 ∼ Bin(n, π 1 ) and y 2 ∼ Bin(n, π 2 ).
- Estimates of the probability of death are given by ˆπ 1 = 1/20 = .05 and ˆπ 2 = 4/20 = .2.
- There are several ways to express the differences in the doses. The advantages and disadvantages of each of the following: 1. Absolute Risk Reduction 2. Relative Risk 3. Odds Ratio 4. Log Odds
Absolute Risk Reduction - this is calculated by taking the difference in the probabilities. The reduction in risk when changing from dose 2 to dose 1 can be estimated by
πˆ 1 − πˆ 2 =. 05 − .2 = −. 15.
The downside is that this difference is difficult to interpret in terms of probability. Note: the difficulty has nothing to do with the negative probability. The ARR will change as we go from dose to dose. This change may or may not be linear - which could make it difficult to describe in terms of covariates. That is, “How is the probability changing from dose to dose as a function (linear or nonlinear) of covariates?”
Relative Risk - this is calculated as a ratio of probabilities. The relative risk of death from dose 1 compared to dose 2 is estimated as
ˆπ 1 /πˆ 2 =. 05 /.2 =. 25.
This is a little more easier to interpret in terms of probabilities. The risk of dying from dose 1 is only 1/4 of the risk of dying from dose 2. Relative risk of less than 1 implies that dose 1 is less poisonous than dose 2. On the downside, a RR of .25 does not say anything about the risk of dying from dose 2. For example, if the probability of dying from dose 1 is .0005 and from dose 2 is .002 then there would be virtually no difference between the doses in practice
Odds Ratio -This measure is calculated as the ratio of odds. An estimate of the ratio of odds is calculated as follows odds(ˆπ 1 ) odds(ˆπ 2 )
The odds of dying from dose 1 is 1/5 versus dose 2, alternatively, dose 2 is 5 times more lethal than dose 1. The downside is the same as for relative risk. An odds ratio of 1 indicates no difference in effects.
Log Odds - This is calculated as the logarithm of the odds ratio. For the budworm example we have as an estimate
log
( odds(ˆπ 1 ) odds(ˆπ 2 )
) = log(.21) = − 1. 56
Note that the log odds varies between −∞ and ∞. A log odds of 0 (i.e. odds ratio of 1) indicates that there are no differences in effect.
Note: The Log Odds measure can also be viewed as the difference of the logits as follows:
log
( odds(π 1 ) odds(π 2 )
) = log
( π 1 /(1 − π 1 ) π 2 /(1 − π 2 )
)
= log
( π 1 1 − π 1
) − log
( π 2 1 − π 2
)
= logit(π 1 ) − logit(π 2 ).
In logistic regression, the difference in logits will correspond to parameters in the logistic regression model.
The Plots (by Gender)
- Proportion Dead vs Dose.
- Empirical logit vs Dose.
- Proportion vs ln(Dose).
- Empirical logit vs ln(Dose).
Note:
- Two candidate models appear reasonable.
- Model 1: π = β 0 + β 1 ln(dose) + β 2 gender + .
- Model 2: logit(π) = β 0 + β 1 gender + β 2 ln(Dose).
- There is no error term on the logit model.
- Probabilities are restricted to the interval [0,1].
- As we will see, Model 1 will have problems.
Model 1 - Simple Linear Regression Approach
- Initial model π = β 0 + β 1 ln(dose) + .
- r^2 for this model is quite high - indicative of a good fit?!
- Predicted values fall within the range of a probability - this is good - but will not always be the case.
- Consider Residuals vs Predicted plot which identifies gender.
- Clearly, gender belongs in the model.
- Reformulated model π = β 0 + β 1 ln(dose) + β 2 gender + .
- r^2 is now at a respectable .98 - indicative of a good fit?!
- Residual patterns look pretty good for this model.
- Notice that one predicted value falls outside the range of [0,1] (Female, Dose=1).
- We will have to abandon this model.
Model 2 - Logistic Regression Approach
The logistic regression model is specified in two steps:
- Random part - specifies the distribution of the observations (not the errors!).
y ∼ Bin(n, π)
- Systematic part - links a function of the expected value of y to a linear combination of the predictors. logit(π) = β 0 + β 1 ln(dose) + β 2 gender
Note:
- The logit function is applied to π itself and not the observed frequencies ˆπ.
- Contrast this to minimizing the least squares in simple linear regression between the response ˆπ and the predictor ln(dose) - where we did apply the observed values to the regression function.
- Technical note: The observed frequencies are applied to the likelihood function which is maximized to obtain parameter estimates.
2.2 Parameter Estimation
Method of Estimation:
There are no closed form solutions to the parameter estimates. The estimates of the parame- ters in logistic regression are obtained by maximizing the likelihood - an iterative procedure. The standard errors, confidence intervals, and hypothesis tests are based on asymptotics (large sample theory). Thus, statistical inference in logistic regression is based on approx- imations. More data is better than less. Technical note: this binomial likelihood has the form of an exponential family which links a function of the mean - nπ using logit as a linear function of regressors.
SAS Implementation Notes: see moths logistic part1.sas
- Descending option - We are modeling the log odds of death as a linear function of logdose and gender2. Without the descending option, SAS will model the probabil- ity of survival (i.e. logit(1 − π)). Why this procedure works this way still alludes my understanding.
- model y/n - y represents the count of dead moths for a particular dose/gender combina- tion, while n represents the total number of moths subjected to a particular dose/gender combination. n must be a variable in the referenced dataset - that is, dividing y by 20 would have generated an error message in the log file.
- gender2 - Moth gender has been converted to 0 for Female and 1 for Male. Converting a dichotomous qualitative variable to 0 or 1 treats gender2 as an indicator variable. As you will see, this will make parameter interpretation easier.
- Suppose we increase ln(dose) by one unit (say from d to d + 1).
- logit(πd+1) = β 0 + β 1 (d + 1) + β 2 (gender) = logit(π) + β 1
- So in terms of the logit(π), β 1 gives the change in log-odds.
- Exponentiating log-odds we get that
odds(πd+1) = odds(π)eβ^1.
- Thus, increasing ln(dose) by one unit increases the odds of death by e^1.^54 = 4.66.
- Note
- The increase in odds is irrespective of the dose from which one increases.
- The increase in the logit is irrespective of the dose from which one increases.
- The increase in odds and in the logit does not depend on gender.
- The probability of death is a non-linear function of d.
- Thus, the probability of death in changing ln(dose) from d to d + 1 does depend on d.
- Compare the fitted logits and probabilities for male moths against ln(dose) in the SAS output.
- The fitted logits are linear in ln(dose).
- The fitted probabilities display an S-curve in ln(dose).
Interpretation of a qualitative predictor.
In the moth data model given by
logit(π) = β 0 + β 1 ln(dose) + β 2 gender,
an estimate for β 2 is given by βˆ 2 ≈ 1 .1.
This parameter can be interpreted in terms of the difference in the logits for males and females. That is,
logit(πF ) − logit(πM ) = β 2 (M ale) − β 2 (F emale) = β 2.
An estimate of β 2 is given by βˆ 2 ≈ 1 .1.
Accordingly, the odds-ratio of dying as a male moth relative to a female moth is given as
e^1.^1 ≈ 3.
This means that the odds of dying for males is 3 times that of females.
Note:
- Direct interpretation of these parameters is possible because of our additive model - that is, we do not have interaction between gender and dose in the model.
- The additive model is based on the assumption that by changing a predictor of interest we can hold the other(s) constant.
Interpretation in the presence of interaction.
Simple interpretations do not hold when the model contains interaction terms of predictors.
Consider adding an interaction term to the moth data. The model can be written as
logit(π) = β 0 + β 1 ln(dose) + β 2 gender + β 3 ln(dose) ∗ gender.
Fitting this model yielded the following estimates.
Analysis of Maximum Likelihood Estimates
Parameter DF Estimate StdError WaldChi-Square Pr > ChiSq Intercept 1 -2.9935 0.5527 29.3353 <. logdose 1 1.3071 0.2411 29.3987 <. gender2 1 0.1752 0.7783 0.0507 0. logdose*gender2 1 0.509 0.3895 1.7078 0.
Analysis of ln(dose) for Females
- Suppose we increase ln(dose) by one unit (say from d to d + 1). We can analyze the difference in the logits for d and d + 1 for females (only) as
- logit(πF, d+1) − logit(πF, d) = β 0 + β 1 (d + 1) − (β 0 + β 1 (d)) = β 1.
- An estimate of β 1 is given by βˆ 1 = 1.3.
- Thus, the odds of death for females is increased by e^1.^3 ≈ 3 .7 for a 1 unit change in ln(dose).
Analysis of ln(dose) for Males
- Suppose we increase ln(dose) by one unit (say from d to d + 1). We can analyze the difference in the logits for d and d + 1 for males (only) as
- logit(πM, d+1)−logit(πM, d) = β 0 +β 1 (d+1)+β 2 +β 3 (d+1)−(β 0 +β 1 (d)+β 2 +β 3 (d)) = β 1 + β 3.
- An estimate of β 1 is given by βˆ 1 = 1.3.
