Chapter 15: A Bernoulli Response

Section A.14 Chapter 15: A Bernoulli Response

Subsection A.14.1 Exercise 15.1

First we save the file “diet.csv" in the working directory of R and read it’s content. Then we apply the function “table" to the two variables in the file in order to obtain the requested frequency table:
🔗
```
diet <- read.csv("_data/diet.csv")
table(diet$health,diet$type)
```
```
##          
##           aha med
##   cancer   15   7
##   death    24  14
##   healthy 239 273
##   illness  25   8
```
The resulting table has two columns and 4 rows. The third row corresponds to healthy subjects. Of these, 239 subjects used the AHA recommended diet and 273 used the Mediterranean diet. We may also plot this data using a mosaic plot:
🔗
```
plot(health~type,data=diet)
```
Examining this plot one may appreciate the fact that the vast majority of the subjects were healthy and the relative proportion of healthy subjects among users of the Mediterranean diet is higher than the relative proportion among users of the AHA recommended diet.
🔗

🔗
In order to test the hypothesis that the probability of keeping healthy following an heart attack is the same for those that use the Mediterranean diet and for those that use the diet recommended by the AHA we create a $2\times 2\text{.}$ This table compares the response of being healthy or not to the type of diet as an explanatory variable. A sequence with logical components, “TRUE" for healthy and “FALSE" for not, is used as the response. Such a sequence is produced via the expression “diet$health==healthy". The table may serve as input to the function “prop.test":
🔗
```
table(diet$health=="healthy",diet$type)
prop.test(table(diet$health=="healthy",diet$type))
```
```
##        
##         aha med
##   FALSE  64  29
##   TRUE  239 273
```
```
## 
##  2-sample test for equality of proportions with continuity
##  correction
## 
## data:  table(diet$health == "healthy", diet$type)
## X-squared = 14.555, df = 1, p-value = 0.0001361
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.1114300 0.3313203
## sample estimates:
##    prop 1    prop 2 
## 0.6881720 0.4667969
```
The function “prop.test" conducts the test that compares between the probabilities of keeping healthy. In particular, the computed $p$-value for the test is $0.0001361\text{,}$ which is less than 0.05. Therefore, we reject the null hypothesis that both diets have the same effect on the chances of remaining healthy following an heart attack.
🔗

🔗
The confidence interval for the difference in probabilities is equal to $[0.1114300, 0.3313203]\text{.}$ The point estimation of the difference between the probabilities is $\hat p_a - \hat p_b = 0.6881720- 0.4667969 \approx 0.22$ in favor of a Mediterranean diet. The confidence interval proposes that a difference as low as 0.11 or as high as 0.33 are not excluded by the data.
🔗

🔗

🔗

Subsection A.14.2 Exercise 15.2

We save the data of the file in a data frame by the name “cushings", produce a mosaic plot with the function “plot" and an histogram with the function “hist":
🔗
```
cushings <- read.csv("_data/cushings.csv")
plot(type~tetra,data=cushings)
hist(cushings$tetra)
```
The mosaic plot describes the distribution of the 4 levels of the response within the different intervals of values of the explanatory variable. The intervals coincide with the intervals that are used in the construction of the histogram. In particular, the third vertical rectangle from the left in the mosaic is associated with the third interval from the left in the histogram[^6]. This interval is associated with the range of values between 20 and 30. The height of the given interval in the histogram is 2, which is the number of patients with “terta" levels that belong to the interval.
🔗

There are 4 shades of *grey* in the first vertical rectangle from the left. Each shade is associated with a different level of the response. The lightest shade of grey, the upmost one, is associated with the level “u". Notice that this is also the shade of grey of the entire third vertical rectangle from the left. The conclusion is that the 2 patients that are associated with this rectangle have Tetrahydrocortisone levels between 2 and 30 and have an unknown type of syndrome.
🔗

🔗
We fit the logistic regression to the entire data in the data frame “cushings" using the function “glm", with the “family=binomial" argument. The response is the indicator that the type is equal to “b". The fitted model is saved in an object called “cushings.fit.all". The application of the function “summary" to the fitted model produces a report that includes the test of the hypothesis of interest:
🔗
```
cushings.fit.all <- glm((type=="b")~tetra,family=binomial, data=cushings)
summary(cushings.fit.all)
```
```
## 
## Call:
## glm(formula = (type == "b") ~ tetra, family = binomial, data = cushings)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0924  -1.0460  -0.8652   1.3426   1.5182  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.12304    0.61330  -0.201    0.841
## tetra       -0.04220    0.05213  -0.809    0.418
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 35.594  on 26  degrees of freedom
## Residual deviance: 34.739  on 25  degrees of freedom
## AIC: 38.739
## 
## Number of Fisher Scoring iterations: 4
```
The test of interest examines the coefficient that is associated wit the explanatory variable “tetra". The estimated value of this parameter is $-0.04220\text{.}$ The $p$-value for testing that the coefficient is 0 is equal to $0.418\text{.}$ Consequently, since the $p$-value is larger than 0.05, we do not reject the null hypothesis that states that the response and the explanatory variable are statistically unrelated.
🔗

Confidence intervals may be computed by applying the function “confint" to the fitted model:
🔗
```
confint(cushings.fit.all)
```
```
##                  2.5 %     97.5 %
## (Intercept) -1.2955624 1.18118256
## tetra       -0.1776113 0.04016772
```
Specifically, the confidence interval for the coefficient that is associated with the explanatory variable is equal to $[-0.1776113, 0.04016772]$
🔗

🔗
If we want to fit the logistic model to a partial subset of the data, say all the observations with values of the response other that “u", we may apply the argument “subset"[^7]. Specifically, adding the expression “subset=(type!=u)" would do the job[^8]. We repeat the same analysis as before. The only difference is the addition of the given expression to the function that fits the model to the data. The fitted model is saved in an object we call “cushings.fit.known":
🔗
```
cushings.fit.known <- glm((type=="b")~tetra,family=binomial, data=cushings,subset=(type!="u"))
summary(cushings.fit.known)
```
```
## 
## Call:
## glm(formula = (type == "b") ~ tetra, family = binomial, data = cushings, 
##     subset = (type != "u"))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2079  -1.1865  -0.7548   1.2033   1.2791  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.11457    0.59947   0.191    0.848
## tetra       -0.02276    0.04586  -0.496    0.620
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 29.065  on 20  degrees of freedom
## Residual deviance: 28.789  on 19  degrees of freedom
## AIC: 32.789
## 
## Number of Fisher Scoring iterations: 4
```
The estimated value of the coefficient when considering only subject with a known type of the syndrome is slightly changed to $-0.02276\text{.}$ The new $p$-value, which is equal to $0.620\text{,}$ is larger than 0.05. Hence, yet again, we do not reject the null hypothesis.
🔗
```
confint(cushings.fit.known)
```
```
##                  2.5 %     97.5 %
## (Intercept) -1.0519135 1.40515473
## tetra       -0.1537617 0.06279923
```
For the modified confidence interval we apply the function “confint". We get now $[-0.1537617, 0.06279923]$ as a confidence interval for the coefficient of the explanatory variable.
🔗

We started with the fitting the model to all the observations. Here we use only the observations for which the type of the syndrome is known. The practical implication of using all observations in the fit is equivalent to announcing that the type of the syndrome for observations of an unknown type is not type “b". This is not appropriate and may introduce bias, since the type may well be “b". It is more appropriate to treat the observations associated with the level “u" as missing observations and to delete them from the analysis. This approach is the approach that was used in the second analysis.
🔗

🔗

🔗

Prev Top Next