STAT 3220 Unit 4 Guided Notes
Unit 4.1-4.2 (Logistic Regression)
Holst, Kromhout, and Brand conducted a case-control study to determine if keeping birds increased the chances of getting lung cancer. It is believed that people who keep birds may inhale more allergens and dust particles and thus may be more likely to get lung cancer. They collected data on 49 people with lung cancer and 98 people with similar demographics who did not have lung cancer. Read more about this study here:
- Gender: F, M: 1 represents females
- Status: High, Low: 1 represents high economic status
- Age: Owner’s age
- Smoked: number of years the person has smoked
- Cigarettes: number of cigarettes smoked per day
- Bird: OwnBird, NoBird: 1 represents own a bird
- Lung Cancer: 1 represents having lung cancer (response)
Step 1: Collect the Data
- Data Compilation: These data were collected from researchers in a reputable journal.
- Data Organization: Are the data organized for logistic regression?- each row contains information about one observation (experimental unit) including all of the explanatory predictors and the response in the columns.
- Suitability of response variable for LOGISTIC regression: We want to look at the breakdown of successes and failures. It is useful if there is a relatively even number of successes and failures in the data as that will allow for more reliable estimates of the coefficients. (Which lead to more accurate predictions.) We do not have formally strict guidelines for this, but it is recommended to have about 10-20 observations per predictor in your model for the lower prevalence group. The bigger discrepancy in ratio of successes to failures, the larger the sample size is needed. We call this the Event Per Variable criterion (EPV).
Defining the response
\(\pi=P(Y=1)=\) The probability of having lung cancer given the set of predictors.
The qualitative variables were originally stored as dummy variables (1 and 0). We are going to rename them as factors and qualitative levels in the data table so that they are easier to interpret in our analysis and output.
# First restore the variable as a factor, instead of numeric
# Then redefine the levels with the ordering of base level then the others
# Response
birdkeeping$LungCancer<-factor(birdkeeping$LungCancer)
levels(birdkeeping$LungCancer)<-c("No","Yes")
# Gender
birdkeeping$Gender<-factor(birdkeeping$Gender)
levels(birdkeeping$Gender)<-c("M","F")
# Status
birdkeeping$Status<-factor(birdkeeping$Status)
levels(birdkeeping$Status)<-c("Low","High")
# Bird
birdkeeping$Bird<-factor(birdkeeping$Bird)
levels(birdkeeping$Bird)<-c("NoBird","OwnBird")# table(response) for counts
head(birdkeeping)table(birdkeeping$LungCancer)
No Yes
98 49 barplot(table(birdkeeping$LungCancer),xlab="Lung Cancer", main="Counts of Response")
Although this break down is not even, there is a suitable number of observations in each group. We have 49 observations in the lower case group, so we would likely not want to include more than about 5 predictors in our model.
Step 2: Hypothesize Relationship (Exploratory Data Analysis)
Quantitative Variables
We explore the relationships with quantitative variables with box plots and means then classify the relationships as existent or not.
- Age: Owner’s age
- Smoked: number of years the person has smoked
- Cigarettes: number of cigarettes smoked per day
# Age
boxplot(Age~LungCancer, data=birdkeeping)
tapply(birdkeeping$Age,birdkeeping$LungCancer,summary)$No
Min. 1st Qu. Median Mean 3rd Qu. Max.
37 52 59 57 63 67
$Yes
Min. 1st Qu. Median Mean 3rd Qu. Max.
37.0 52.0 58.0 56.9 63.0 66.0 # Smoked
boxplot(Smoked~LungCancer, data=birdkeeping)
tapply(birdkeeping$Smoked,birdkeeping$LungCancer,summary)$No
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 16.00 25.00 24.96 38.75 47.00
$Yes
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 27.00 36.00 33.63 41.00 50.00 # Cigarettes
boxplot(Cigarettes~LungCancer, data=birdkeeping)
tapply(birdkeeping$Cigarettes,birdkeeping$LungCancer,summary)$No
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 6.00 15.00 14.21 20.00 45.00
$Yes
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 15.00 20.00 18.82 20.00 40.00 - We do not see an association with the age of a person smoked and having lung cancer. For this sample, the average age is about 57 years old for both those who had and did not have lung cancer.
- We see a slight association with the number of years a person smoked and having lung cancer. On average, those who have lung cancer smoked longer.
- We see a slight association with the number of cigarettes a person smoked per day and having lung cancer. On average, those who have lung cancer smoked more cigarettes per day (20, compared to 15)
Qualitative Variables
We will explore the relationship with qualitative variables with mosaic plots and conditional proportions and classify each relationship as existent or not.
- Gender: M, F
- Status: High, Low
- Bird: NoBird, OwnBird
# Gender
#Mosaic Plots and proportions for qualitative variables
mosaicplot(LungCancer~Gender,data=birdkeeping)
#contingency table for counts (table(rows, columns))
table(birdkeeping$Gender,birdkeeping$LungCancer)
No Yes
M 74 37
F 24 12#Convert contingency table to proportions, margin=2 specifies totaling over the columns
prop.table(table(birdkeeping$Gender,birdkeeping$LungCancer),margin=2)
No Yes
M 0.755102 0.755102
F 0.244898 0.244898# Status
#Mosaic Plots and proportions for qualitative variables
mosaicplot(LungCancer~Status,data=birdkeeping)
#contingency table for counts (table(rows, columns))
table(birdkeeping$Status,birdkeeping$LungCancer)
No Yes
Low 65 37
High 33 12#Convert contingency table to proportions, margin=2 specifies totaling over the columns
prop.table(table(birdkeeping$Status,birdkeeping$LungCancer),margin=2)
No Yes
Low 0.6632653 0.7551020
High 0.3367347 0.2448980# Bird
#Mosaic Plots and proportions for qualitative variables
mosaicplot(LungCancer~Bird,data=birdkeeping)
#contingency table for counts (table(rows, columns))
table(birdkeeping$Bird,birdkeeping$LungCancer)
No Yes
NoBird 64 16
OwnBird 34 33#Convert contingency table to proportions, margin=2 specifices totalling over the columns
prop.table(table(birdkeeping$Bird,birdkeeping$LungCancer),margin=2)
No Yes
NoBird 0.6530612 0.3265306
OwnBird 0.3469388 0.6734694We do not see an association between lung cancer and gender. The proportion of males (or the proportion of females) is exactly the same for both those who have and not have lung cancer.
We do not see an association between lung cancer and status. The proportion of high (or the proportion of low) is about the same for both those who have and not have lung cancer.
We see there might be an association with whether a person owned a bird and whether they have lung cancer. We that for those without lung cancer, 35 percent owned a bird. Compared to the 67 percent of those who had lung cancer owning a bird.
We will not explore interactions in this example.
We can also evaluate multicollinearity between the predictors.
We can also use variable screening in logistic regression.
Step 3: Estimate the model parameters
Here we will fit the model with all of the predictors.
NOTE: It is good practice to define what you want to be treated as the base level, or “failure” Y=0. You can do this with the relevel function. The relevel function can be used even if the original column is numeric.
# Be sure to clearly define the base level (failure group) for response
birdkeeping$LungCancer<-relevel(factor(birdkeeping$LungCancer),"No")
# We can use the . to abbreviate the syntax and fit a model with all
# predictors in the data table.
#notice birdmod1 and 1a result in the same model
birdmod1<-glm(LungCancer~ Gender+ Status+ Age+Smoked+Cigarettes+Bird,data=birdkeeping,family=binomial())
summary(birdmod1)
Call:
glm(formula = LungCancer ~ Gender + Status + Age + Smoked + Cigarettes +
Bird, family = binomial(), data = birdkeeping)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.93736 1.80425 -1.074 0.282924
GenderF 0.56127 0.53116 1.057 0.290653
StatusHigh 0.10545 0.46885 0.225 0.822050
Age -0.03976 0.03548 -1.120 0.262503
Smoked 0.07287 0.02649 2.751 0.005940 **
Cigarettes 0.02602 0.02552 1.019 0.308055
BirdOwnBird 1.36259 0.41128 3.313 0.000923 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 187.14 on 146 degrees of freedom
Residual deviance: 154.20 on 140 degrees of freedom
AIC: 168.2
Number of Fisher Scoring iterations: 5birdmod1a<-glm(LungCancer~.,data=birdkeeping,family=binomial())
summary(birdmod1a)
Call:
glm(formula = LungCancer ~ ., family = binomial(), data = birdkeeping)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.93736 1.80425 -1.074 0.282924
GenderF 0.56127 0.53116 1.057 0.290653
StatusHigh 0.10545 0.46885 0.225 0.822050
Age -0.03976 0.03548 -1.120 0.262503
Smoked 0.07287 0.02649 2.751 0.005940 **
Cigarettes 0.02602 0.02552 1.019 0.308055
BirdOwnBird 1.36259 0.41128 3.313 0.000923 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 187.14 on 146 degrees of freedom
Residual deviance: 154.20 on 140 degrees of freedom
AIC: 168.2
Number of Fisher Scoring iterations: 5MODEL:
- Log-Odds: \(\pi^*=\ln(\frac{\pi}{1-\pi})=\beta_0+\beta_1 GenderF+\beta_2 StatusHigh+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6BirdOwnBird\)
- Odds: \(\frac{\pi}{1-\pi}=e^{\beta_0+\beta_1 GenderF+\beta_2 StatusHigh+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6BirdOwnBird}\)
- Probability of Success: \(\pi=\dfrac{e^{\beta_0+\beta_1 GenderF+\beta_2 StatusHigh+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6BirdOwnBird}}{1+ e^{\beta_0+\beta_1 GenderF+\beta_2 StatusHigh+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6BirdOwnBird}}\)
PREDICTION EQUATION:
- Estimated Log-Odds:\(\hat{\pi^*}=-1.937+0.56 GenderF+0.105 StatusHigh-0.039 Age+0.073 Smoked+0.026Cigarettes+1.36BirdOwnBird\)
- Estimated Odds:\(\widehat{\frac{\pi}{1-\pi}}=e^{-1.937+0.56 GenderF+0.105 StatusHigh-0.039 Age+0.073 Smoked+0.026Cigarettes+1.36BirdOwnBird}\)
- Predicted Probability of Success:\(\hat{\pi}=\dfrac{e^{-1.937+0.56 GenderF+0.105 StatusHigh-0.039 Age+0.073 Smoked+0.026Cigarettes+1.36BirdOwnBird}}{1+e^{-1.937+0.56 GenderF+0.105 StatusHigh-0.039 Age+0.073 Smoked+0.026Cigarettes+1.36BirdOwnBird}}\)
Where:
- GenderF= 1 if female, 0 otherwise
- StatusHigh= 1 if high economic status, 0 otherwise
- BirdOwnBird= 1 if own a bird, 0 otherwise
Math Interlude
Show that the beta coefficients do not change with whatever form of the model you write.
Consider the model \(\pi^*=ln(\frac{\pi}{1-\pi})=\beta_0+\beta_1x\)
Solve the equation for \(\pi\) to show \(\pi=\dfrac{e^{\beta_0+\beta_1x}}{1+e^{\beta_0+\beta_1x}}\). Show all of your steps. Use the appropriate markdown syntax (it may help to to write it out on paper before writing into the markdown file.)
Show a one unit change in x results in a \(e^\beta\) factor increase in the odds for quantitative variables.
First, consider the model for the log-odds: \(\pi^*=\beta_0+\beta_1x\). Compute \(\pi^*\) for x=0, x=1, and x=2 in terms of the betas. Then determine how much do the log-odds change for a one unit change in x.
| X | Log-odds: \(\pi^*=\beta_0+\beta_1x\) |
|---|---|
| X=0 | \(=\beta_0\) |
| X=1 | \(=\beta_0+\beta_1\) |
| X=2 | \(=\beta_0+2*\beta_1\) |
For each additional unit increase in x, the log-odds increases by \(\beta_1\). (Similiar to how we interpret slope in MLR)
Now, consider the model for the odds:\(\dfrac{\pi}{1-\pi}=e^{\beta_0+\beta_1x}\). Compute the odds for x=0, x=1, and x=2 in terms of the betas. Then determine how much do the odds change for a one unit change in x.
| X | odds: \(\frac{\pi}{1-\pi}=e^{\beta_0+\beta_1x}\) |
|---|---|
| X=0 | \(=e^{\beta_0}\) |
| X=1 | \(=e^{\beta_0+\beta_1}=e^{\beta_0}*e^{\beta_1}\) |
| X=2 | \(=e^{\beta_0+2*\beta_1}=e^{\beta_0}*e^{2\beta_1}=e^{\beta_0}*e^{\beta_1}*e^{\beta_1}\) |
For each additional unit increase in x, the odds increases by a multiplicative factor of \(e^{\beta_1}\).
Show the interpretation for a dummy variable coefficient.
Assume we have a qualitative variable with 3 levels (A, B, C). Let DummyA = 1 if group A and 0 otherwise and DummyB=1 if group B, 0 otherwise. Group C is the base level.
Consider the model for the odds:\(\dfrac{\pi}{1-\pi}=e^{\beta_0+\beta_1DummyA +\beta_2DummyB}\). Compute the odds group A, B, and C in terms of the betas. Then determine the interpretation of a single beta coefficient associated with a dummy variable.
| Group | odds: \(\frac{\pi}{1-\pi}=e^{\beta_0+\beta_1DummyA+\beta_2DummyB}\) |
|---|---|
| Odds of Group A | \(=e^{\beta_0+\beta_1}\) |
| Odds of Group B | \(=e^{\beta_0+\beta_2}\) |
| Odds of Group C | \(=e^{\beta_0}\) |
| Coefficient | odds |
|---|---|
| \(=e^{\beta_0}\) | Odds of Group C (the base level) |
| \(=e^{\beta_1}\) | Odds of Group A / Odds of Group C |
| \(=e^{\beta_2}\) | Odds of Group B/ Odds of Group C |
INTERPRETATIONS
# Extract just coefficients (or intervals)
coefficients(birdmod1)(Intercept) GenderF StatusHigh Age Smoked Cigarettes
-1.93735861 0.56127270 0.10544761 -0.03975542 0.07286848 0.02601689
BirdOwnBird
1.36259456 confint(birdmod1, level=.95)Waiting for profiling to be done... 2.5 % 97.5 %
(Intercept) -5.59554204 1.54015158
GenderF -0.48194383 1.61625027
StatusHigh -0.82695901 1.02585600
Age -0.11263170 0.02827386
Smoked 0.02529552 0.13082279
Cigarettes -0.02407110 0.07680737
BirdOwnBird 0.57309984 2.19296212#Exponentiation of betas (or intervals)
exp(coefficients(birdmod1))(Intercept) GenderF StatusHigh Age Smoked Cigarettes
0.1440840 1.7529020 1.1112079 0.9610245 1.0755891 1.0263583
BirdOwnBird
3.9063153 exp(confint(birdmod1,level=.95))Waiting for profiling to be done... 2.5 % 97.5 %
(Intercept) 0.003714385 4.665297
GenderF 0.617581748 5.034178
StatusHigh 0.437377326 2.789482
Age 0.893479671 1.028677
Smoked 1.025618164 1.139766
Cigarettes 0.976216301 1.079834
BirdOwnBird 1.773756897 8.961720Quantitative Interpretations:
- \(\hat{\beta_3}\): For a one unit year in age, the estimated log-odds (\(\hat{\pi^*}\)) of a person having lung cancer decreases by 0.034, given the other variables are held constant.
- \(e^{\hat{\beta_3}}\): For a one unit year in age, the estimated odds (\(\widehat{\frac{\pi}{1-\pi}}\)) of a person having lung cancer decreases by a factor of 0.96, given the other variables are held constant.
- Confidence interval for \(e^{\beta_3}\): We are 95% confident that for a one unit year in age, the odds of a person having lung cancer changes by a factor between 0.893 and 1.02, given the other variables are held constant. (Since 1 is in this interval the variable “age” is not a significant predictor.)
On Your Own: Try the interpretations for the “Smoked” variable.
Qualitative Interpretations:
- \(\hat{\beta_6}\): The estimated log-odds (\(\hat{\pi^*}\)) of a person having lung cancer is 1.36 higher for those who own a bird, given the other variables are held constant.
- \(e^{\hat{\beta_6}}\): The estimated odds (\(\widehat{\frac{\pi}{1-\pi}}\)) of a person having lung cancer is 3.9 times higher for those who own a bird, given the other variables are held constant.
- Confidence interval for \(e^{\beta_3}\): We are 95% confident that the odds of a person having lung cancer increases by a factor between 1.77 and 8.96 given the other variables are held constant. (Since 1 is not in this interval “bird” is a significant predictor)
On Your Own: Try the interpretations for the “Status” variable.
Notes on Intervals and Significance
For logistic regression we use the hypothesis format of \(H_0: \beta_i=0\). That translates to \(H_0: e^{\beta_i}=1\). Therefore, if and only if the interval of the odds does not contain 1, the interval for the log-odds will not contain 0, and the parameter is significant.
Therefore, based on the intervals, we can see that Smoked and Bird are significant. Their intervals for log-odds do not contain 0 and their intervals for odds do not contain 1.
Step 4: Evaluate the model
\(\pi^*=\ln(\frac{\pi}{1-\pi})=\beta_0+\beta_1 Gender+\beta_2 Status+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6Bird\)
First we Perform the model adequacy test (Likelihood Ratio Test):
#REPRINT MODEL
birdmod1<-glm(LungCancer~.,data=birdkeeping,family=binomial())
summary(birdmod1)
Call:
glm(formula = LungCancer ~ ., family = binomial(), data = birdkeeping)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.93736 1.80425 -1.074 0.282924
GenderF 0.56127 0.53116 1.057 0.290653
StatusHigh 0.10545 0.46885 0.225 0.822050
Age -0.03976 0.03548 -1.120 0.262503
Smoked 0.07287 0.02649 2.751 0.005940 **
Cigarettes 0.02602 0.02552 1.019 0.308055
BirdOwnBird 1.36259 0.41128 3.313 0.000923 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 187.14 on 146 degrees of freedom
Residual deviance: 154.20 on 140 degrees of freedom
AIC: 168.2
Number of Fisher Scoring iterations: 5#Compute the test statistic
lltest<-birdmod1$null.deviance-birdmod1$deviance
lltest[1] 32.93675#compute the pvalue
#the second argument is the degrees of freedom = k
llpvalue<-1-pchisq(lltest,6)
llpvalue[1] 1.078379e-05- Hypotheses:
- \(H_0: \beta_1=...=\beta_6=0\) (the model is not adequate)
- \(H_a\):at least one of \(\beta_1 ,...,\beta_6\) does not equal 0 (the model is adequate)
- Distribution of test statistic: chi-square with 6 DF
- Test Statistic: \(\chi^2=187.14-154.20=32.93675\)
- Pvalue: <0.0001
- Decision: 0.0001<0.05 -> REJECT H0
- Conclusion: The model with these predictors is adequate at predicting the probability of having lung cancer.
Then we test “the most important predictors”:
Test the Significance of Bird (one predictor- Wald Test)
- Hypotheses:
- \(H_0: \beta_6=0\)
- \(H_a:\beta_6 \neq 0\)
- Distribution of test statistic: Chi Square with 1 DF OR Z
- Test Statistic: \(\chi^2=3.31^2\) OR Z=3.31
- Pvalue: 0.000923
- Decision: 0.000923<0.05 -> REJECT H0
- Conclusion: Owning a bird contributes to predicting the probability of having lung cancer, given the other variables are present in the model. (Additionally, since 1 was not in the interval for the odds “bird” is a significant predictor.)
Test the Significance of Age (one predictor- Wald Test)
- Hypotheses:
- \(H_0: \beta_3=0\)
- \(H_a:\beta_3 \neq 0\)
- Distribution of test statistic: Chi Square with 1 DF OR Z
- Test Statistic: \(\chi^2=|-1.2544|^2\) OR Z=-1.120
- Pvalue: 0.262503
- Decision: 0.262503<0.05 -> FAIL TO REJECT H0
- Conclusion: Age does not contribute to predicting the probability of having lung cancer, given the other variables are present in the model. (Additionally, since 1 was in the interval for the odds “age” is not a significant predictor.) We will remove it and refit the model.
birdmod2<-glm(LungCancer~.-Age,data=birdkeeping,family=binomial())
summary(birdmod2)
Call:
glm(formula = LungCancer ~ . - Age, family = binomial(), data = birdkeeping)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -3.80143 0.82243 -4.622 3.8e-06 ***
GenderF 0.58361 0.52725 1.107 0.268340
StatusHigh 0.03958 0.46451 0.085 0.932102
Smoked 0.05677 0.02033 2.793 0.005229 **
Cigarettes 0.03107 0.02476 1.255 0.209555
BirdOwnBird 1.43560 0.40655 3.531 0.000414 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 187.14 on 146 degrees of freedom
Residual deviance: 155.49 on 141 degrees of freedom
AIC: 167.49
Number of Fisher Scoring iterations: 5Step 5: Assess the model (Unit 4.2)
#Compute AIC and BIC
AIC(birdmod2)[1] 167.4945BIC(birdmod2)[1] 185.4371#Hosmer Lemeshow test
#Requires glmtoolbox package
hltest(birdmod2)
The Hosmer-Lemeshow goodness-of-fit test
Group Size Observed Expected
1 15 0 0.6212571
2 15 3 1.5580862
3 15 2 2.2771044
4 15 2 3.2741449
5 15 6 4.1286036
6 15 5 4.9967048
7 15 5 6.2345576
8 15 8 8.2309190
9 15 9 9.2164032
10 12 9 8.4622194
Statistic = 4.54343
degrees of freedom = 8
p-value = 0.80507 We will evaluate the logistic regression model with AIC, BIC, and Hosmer-Lemeshow Goodness of Fit Test. Individually interpreting these metrics are not practically useful. But these are more typically used to compare between two or more models to determine which is best.
- AIC = 167.49
- BIC = 185.44
- Hosmer-Lemeshow Test:
- \(H_0\): The logistic regression model provides an adequate fit to the data
- With a p-value of 0.805 , we fail to reject the null hypothesis. There is not enough evidence to determine the model does not provide a good fit. We can move forward in the analysis.
Step 6: Check Residuals and Assumptions
We will not cover this specifically in this course, but the supplemental material gives you a bit of a background.
Step 7: Use the model for prediction
coef(birdmod2)(Intercept) GenderF StatusHigh Smoked Cigarettes BirdOwnBird
-3.80143198 0.58361083 0.03957626 0.05677305 0.03106926 1.43559802 - Model: \(\pi=\dfrac{e^{\beta_0+\beta_1 Gender+\beta_2 Status+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6Bird}}{1+e^{\beta_0+\beta_1 Gender+\beta_2 Status+\beta_3 Age+\beta_4 Smoked+\beta_5Cigarettes+\beta_6Bird}}\)
- Prediction Equation: \(\hat{\pi}=\dfrac{e^{-3.801 +0.584GenderF+ 0.0396StatusHigh+ 0.057Smoked +0.031Cigarettes+ 1.436BirdOwnBird}}{1+e^{-3.801 +0.584GenderF+ 0.0396StatusHigh+ 0.057Smoked +0.031Cigarettes+ 1.436BirdOwnBird}}\)
Where:
- GenderF= 1 if female, 0 otherwise
- StatusHigh= 1 if high economic status, 0 otherwise
- BirdOwnBird= 1 if own a bird, 0 otherwise
#compute pi-hat values. These are probabilities of success for each value.
birdkeeping$fits<-fitted.values(birdmod2)
head(birdkeeping)#summarize the "correctly" classified predictions
blr_confusion_matrix(birdmod2)Confusion Matrix and Statistics
Reference
Prediction No Yes
0 81 23
1 17 26
Accuracy : 0.7279
No Information Rate : 0.6667
Kappa : 0.3684
McNemars's Test P-Value : 0.4292
Sensitivity : 0.5306
Specificity : 0.8265
Pos Pred Value : 0.6047
Neg Pred Value : 0.7788
Prevalence : 0.3333
Detection Rate : 0.1769
Detection Prevalence : 0.2925
Balanced Accuracy : 0.6786
Precision : 0.6047
Recall : 0.5306
'Positive' Class : 1- Interpreting \(\hat{\pi}}\): For a Male with low economic class, who smoked for 19 years, 12 cigarettes, owned a bird has a 0.286 probability of having lung cancer.
- We can compare how well our model does at predicting within the data set. For example, observation 1 was actually a “success”, but had a predicted probability of 0.28.
Final note: We can classify “correct” and “incorrect” conclusions. The default cutoff is typically 0.5. Meaning is \(\hat{\pi}\) if greater than 0.5 it will be labeled as a “success”. You can change this threshold depending on the context.
- accuracy overall the percentage of correctly identified results. (0.7279)
- sensitivity is the ability of the model to correctly identify a success (true positive rate) (0.5306)
- specificity is the ability of the model to correctly identify a failure (true negative rate) (0.8265)
- McNemars’s Test has a null that there is no difference
between the proportion of false positives and false negatives (good
thing). A pvalue greater than alpha means our model doesn’t over correct
in either direction.
- Here we have a p-value of 0.4292. So that is good for our assessment. Despite the rates of false positive
4.3 One Factor ANOVA
A sample of 44 healthy male college students participated in an experiment, Each student was asked to memorize a list of 40 words (20 on a green list and 20 on a red list).The students were then randomly assigned of one of 4 treatments: Group A: received 2 alcoholic drinks, Group AC: received 2 alcoholic drinks with caffeine powder dissolved in their drinks, Group AR: received 2 alcoholic drinks and a monetary reward for correct responses on the task, Group P: were told they received 2 alcoholic drinks, but instead just received a carbonated drink.
After consuming their drinks and resting for 25 minutes, the students performed the word completion task. The response variable ”task score” represents the difference between the proportion of correct responses on the green list and incorrect on the red list.
Step 1: Collect the Data and Determine Experimental Design
Consider the design of the data collection: For the One and Two Factor ANOVA we are covering the analysis on completely randomized designs. We want to ensure our data meet that criteria and define the other components of the experimental design.
drinkersdata<-read.delim("https://raw.githubusercontent.com/kvaranyak4/STAT3220/main/DRINKERS.txt",header=T)
names(drinkersdata)[1] "GROUP" "SCORE"head(drinkersdata)table(drinkersdata$GROUP)
A AC AR P
11 11 11 11 Components of Experimental design
This is a completely randomized design because subjects are randomly assigned to treatments. There is one factor, group, with four levels: AR, AC, A, and P. Therefore, there are four treatments: AR, AC, A, P. It is balanced because each treatment has the same number of experimental units. The response variable is the task score.
Step 2: Hypothesize Relationship (Exploratory Data Analysis)
We will look at side by side box plots of the treatment groups and compare the mean value of the response. We can also compare the 5 number summary for a more detailed interpretation. (We could also determine the approximate distribution of each treatment group- in the next section we will see that we want the treatments to have normally distributed data.)
boxplot(SCORE~GROUP, data=drinkersdata)
tapply(drinkersdata$SCORE,drinkersdata$GROUP,summary)$A
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.35000 -0.05000 0.16000 0.06364 0.19000 0.31000
$AC
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0200 0.1750 0.2200 0.2655 0.3750 0.5000
$AR
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 0.400 0.500 0.440 0.525 0.610
$P
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.12 0.33 0.43 0.40 0.47 0.62 - It appears there are differences in the mean task scores for the four groups, so we will perform a statistical analysis to compare the means through the ANOVA procedure.
Step 3: Perform the ANOVA
- Null Hypothesis \(H_0 : \mu_A = \mu_{AC} = \mu_{AR} = \mu_P\)
- Alternative Hypothesis at least two means are different
drinkanova1<-aov(SCORE~GROUP, data=drinkersdata)
summary(drinkanova1) Df Sum Sq Mean Sq F value Pr(>F)
GROUP 3 0.9506 0.3169 10.29 3.76e-05 ***
Residuals 40 1.2317 0.0308
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Distribution of Test Statistic : F with numerator degrees of freedom, v1 = 4 - 1 = 3 and denominator degrees of freedom, v2 = 44 - 4 = 40 (F with 3, 40)
- Test Statistic F test: F = MST/MSE = 0.317/0.0308 = 10.29 with pvalue= 3.76e-05
- CONCLUSION: With a pvalue this small, we reject the null hypothesis. There is evidence that at least two population treatment means are different from each other. Based on the sample means, it appears the group with just alcohol has a much lower than the rest, while the group with a monetary reward had a much higher mean than the rest. We will perform the post-hoc test to compare pairwise means in the next section.
4.3 Two Factor ANOVA
The chemical element antimony is sometimes added to tin – lead solder to replace the more expensive tin and to reduce the cost of soldering. A factorial experiment was conducted to determine how antimony affects the strength of the tin–lead solder joint (Journal of Materials Science, May 1986). Tin–lead solder specimens were prepared using one of four possible cooling methods (water- quenched, WQ; oil-quenched, OQ; air-blown, AB; and furnace-cooled, FC) and with one of four possible amounts of antimony (0%, 3%, 5%, and 10%) added to the composition. Three solder joints were randomly assigned to each of the 4 × 4 = 16 treatments and the shear strength of each measured.
Step 1: Collect the Data and Determine Experimental Design
Consider the design of the data collection: For the One and Two Factor ANOVA we are covering the analysis on completely randomized designs. We want to ensure our data meet that criteria and define the other components of the experimental design.
tindata<-read.delim("https://raw.githubusercontent.com/kvaranyak4/STAT3220/main/TINLEAD.txt",header=T)
names(tindata)[1] "ANTIMONY" "METHOD" "STRENGTH"head(tindata)table(tindata$ANTIMONY,tindata$METHOD)
AB FC OQ WQ
0 3 3 3 3
3 3 3 3 3
5 3 3 3 3
10 3 3 3 3Components of Experimental design
This is a completely randomized design because experimental units are randomly assigned to treatments. There two factors, cooling method, with four levels: AB, FC, OQ, and WQ and antimony with four levels 0%, 3%, 5%, and 10%. Therefore, there are 16 treatments. It is balanced because each treatment has the same number of experimental units. The response variable is the sheer strength.
Step 2: Hypothesize Relationship (Exploratory Data Analysis)
For two factor designs, we will examine the interaction plot first.
# We want to begin with examining the interaction
interaction.plot(tindata$METHOD, tindata$ANTIMONY, tindata$STRENGTH,fun=mean,trace.label="Antimony Level"
, xlab="Cooling Method",ylab="Mean Sheer Strength",main="Interaction Plot Cooling Method X Antimony")
# We can also look at the main effects
#cooling method
boxplot(STRENGTH~METHOD, data=tindata)
tapply(tindata$STRENGTH,tindata$METHOD,summary)$AB
Min. 1st Qu. Median Mean 3rd Qu. Max.
15.80 18.05 20.70 20.18 22.30 22.90
$FC
Min. 1st Qu. Median Mean 3rd Qu. Max.
16.40 17.60 19.50 18.95 20.00 20.90
$OQ
Min. 1st Qu. Median Mean 3rd Qu. Max.
16.40 19.75 20.50 20.45 21.15 24.30
$WQ
Min. 1st Qu. Median Mean 3rd Qu. Max.
15.20 17.48 18.80 18.64 19.50 22.30 #antimony level
boxplot(STRENGTH~ANTIMONY, data=tindata)
tapply(tindata$STRENGTH,tindata$ANTIMONY,summary)$`0`
Min. 1st Qu. Median Mean 3rd Qu. Max.
17.60 19.12 19.80 20.18 20.70 24.30
$`3`
Min. 1st Qu. Median Mean 3rd Qu. Max.
18.60 19.38 20.20 20.41 21.10 22.90
$`5`
Min. 1st Qu. Median Mean 3rd Qu. Max.
16.40 19.68 20.55 20.62 21.77 22.90
$`10`
Min. 1st Qu. Median Mean 3rd Qu. Max.
15.20 16.40 17.10 17.02 17.60 19.00 - INTERACTION: While there is some crossing in the interaction plot, there is not strong evidence of an interaction because generally the hierarchy stays the same across the levels. Additionally there are very small sample sizes in each treatment (3), so likely an interaction will not be significant.
- MAIN EFFECT (METHOD): The average value of the response appears similar across all groups.
- MAIN EFFECT (ANTIMONY): The average value of the response appears similar for 0%, 3%, and 5%, but different for 10%.
Step 3: Perform the ANOVA
Significance of Interaction
NOTE: With a two factor factorial design, we will always test the significance of the interaction first.
- Null Hypothesis \(H_0 : \mu_{AB,0} = \mu_{FC,0} = \mu_{OQ,0} = \mu_{WQ,0}=...=\mu_{WQ,10}\) (there are 16 means!)
- Alternative Hypothesis at least two means are different
# We include the interaction
# Because Antimony was read in as a numeric variable, we want to
# indicate that it is actually a factor
tinanova1<-aov(STRENGTH~METHOD*factor(ANTIMONY), data=tindata)
summary(tinanova1) Df Sum Sq Mean Sq F value Pr(>F)
METHOD 3 28.63 9.54 5.527 0.00357 **
factor(ANTIMONY) 3 104.19 34.73 20.117 1.64e-07 ***
METHOD:factor(ANTIMONY) 9 25.13 2.79 1.617 0.15226
Residuals 32 55.25 1.73
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Distribution of Test Statistic: F with (4-1)(4-1), 4*4(3-1)= 9, 32 DF
- Test Statistic F test: F = 1.62 with pvalue 0.15226
- CONCLUSION: With a pvalue this large, we fail to reject the null hypothesis. There is not evidence of an interaction. This means that the shear strength for each level of cooling method does not depend on the antimony level.
- Next step: remove the interaction and test the main effect
Significance of Main effects
#We do not include the interaction, just add the terms
# Because Antimony was read in as a numeric variable, we want to
# indicate that it is actually a factor
tinanova2<-aov(STRENGTH~METHOD+factor(ANTIMONY), data=tindata)
summary(tinanova2) Df Sum Sq Mean Sq F value Pr(>F)
METHOD 3 28.63 9.54 4.868 0.00549 **
factor(ANTIMONY) 3 104.19 34.73 17.716 1.58e-07 ***
Residuals 41 80.38 1.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Test for Method:
- Null Hypothesis \(H_0 : \mu_{AB} = \mu_{FC} = \mu_{OQ} = \mu_{WQ}\)
- Alternative Hypothesis at least two means are different
- Test Statistic: 4.87 with p-value 0.0055
- Conclusion: We can conclude at least one of the types of cooling methods produces a different shear strengths of the solders than the others.
- Test for Antimony:
- Null Hypothesis \(H_0 : \mu_0 = \mu_3 = \mu_5 = \mu_{10}\)
- Alternative Hypothesis at least two means are different
- Test Statistic: 17.72 with p-value <.0001
- Conclusion: We can conclude at least one of the levels of antimony produces a different shear strengths of the solders than the others.
4.4 POST HOC + ASSUMPTIONS
Drinkers Example Continued
Step 4 Post Hoc Analysis
The question of interest is: Does coffee or some other form of stimulation really allow a person suffering from alcohol intoxication to ”sober up”? What is your conclusion? Include any relevant inferences and confidence intervals. Use alpha=0.05.
TukeyHSD(drinkanova1,conf.level=.95) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = SCORE ~ GROUP, data = drinkersdata)
$GROUP
diff lwr upr p adj
AC-A 0.2018182 0.001256169 0.4023802 0.0480809
AR-A 0.3763636 0.175801623 0.5769256 0.0000618
P-A 0.3363636 0.135801623 0.5369256 0.0003282
AR-AC 0.1745455 -0.026016558 0.3751075 0.1075832
P-AC 0.1345455 -0.066016558 0.3351075 0.2892292
P-AR -0.0400000 -0.240562013 0.1605620 0.9501001plot(TukeyHSD(drinkanova1,conf.level=.95))
# Compute Minimum Significant Difference w for interactions
# q(p,v) Critical studentized range
q<-qtukey(.95,nmeans=4,df=40)
q[1] 3.790685s<-sqrt(0.0308)
n<-11
w<-q*(s/sqrt(n))
w[1] 0.2005842IF AND ONLY IF - The absolute value of the difference (diff) is greater than w - Then, the confidence interval does not contain 0 - Then the pvalue for Tukey’s test is less than alpha, - Then the pair of means are significantly different.
We see for Caffeine- Alcohol (AC-A) - The |difference|= |0.2018182| is greater than w=0.2005842 - The interval (0.001256169, 0.4023802) does not contain 0. - The pvalue = 0.0480809 is less than alpha=0.05 - Therefore the mean difference in recall scores between the group with caffeine is significantly different from the group with just Alcohol.
We can also look at the confidence intervals for the difference in the mean scores of the alcohol group compared to the others. Note all of the confidence intervals are positive. This implies that the alcohol group under performs compared to every other group.
| Pair | Lower | Upper |
|---|---|---|
| AC-A | 0.00126 | 0.4024 |
| AR-A | 0.1758 | 0.5769 |
| P-A | 0.1358 | 0.5369 |
Interpretation of Confidence Interval: We are 95% confident that the true mean difference of scores for those who also had caffeine is between 0.00126 and 0.40 higher than the group with just alcohol. There is a significant difference between the group with just alcohol and the group with caffeine, but it may not be practically significant.
Step 5: Assumptions
Normality
#NOTE: this function is qqPlot from the car package
qqPlot(drinkersdata$SCORE, groups=drinkersdata$GROUP)
We see the normality assumption is reasonably met because the scores of each treatment group follow closely on the line of the QQ Plot. We will use Bartlett’s Test of Homogeneity.
Equal Variance
Because we determined the normality assumption was reasonably met, we will use Bartlett’s test of homogeneity.
#this function is in the stats package
bartlett.test(SCORE~GROUP, data=drinkersdata)
Bartlett test of homogeneity of variances
data: SCORE by GROUP
Bartlett's K-squared = 1.7548, df = 3, p-value = 0.6248#Visualize Constant Variance
boxplot(SCORE~GROUP, data=drinkersdata)
- The null hypothesis is that the variance of each treatment is equal. With a p-value this large, we determine the equal variance assumption is met for these data effects.
Tin Example Continued
Step 4: Post Hoc Analysis (multiple comparisons test)
Since the main effects are significant, we go to post hoc analysis (Tukey’s test in this class).
TukeyHSD(tinanova2,conf.level = 0.95) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = STRENGTH ~ METHOD + factor(ANTIMONY), data = tindata)
$METHOD
diff lwr upr p adj
FC-AB -1.2250000 -2.75555373 0.305553731 0.1567570
OQ-AB 0.2750000 -1.25555373 1.805553731 0.9628573
WQ-AB -1.5333333 -3.06388706 -0.002779602 0.0494355
OQ-FC 1.5000000 -0.03055373 3.030553731 0.0565864
WQ-FC -0.3083333 -1.83888706 1.222220398 0.9488521
WQ-OQ -1.8083333 -3.33888706 -0.277779602 0.0149509
$`factor(ANTIMONY)`
diff lwr upr p adj
3-0 0.2333333 -1.297220 1.763887 0.9767248
5-0 0.4416667 -1.088887 1.972220 0.8663147
10-0 -3.1583333 -4.688887 -1.627780 0.0000119
5-3 0.2083333 -1.322220 1.738887 0.9832100
10-3 -3.3916667 -4.922220 -1.861113 0.0000031
10-5 -3.6000000 -5.130554 -2.069446 0.0000010plot(TukeyHSD(tinanova2,conf.level = 0.95))
- Method:
- Water Quenched and Oil Quenched are the only pairwise means that are significantly different from each other. (WQ-AB have a pvalue very close to alpha- what does that mean?? We need to interpret pvalues!)
- We are 95% confident that the Water Quenched solders have a shear strength between 0.277 and 3.339 lower than the oil quenched group.
- Antimony:
- The 10% group is significantly different from the rest, but the other groups are not different from each other.
- We are 95% confident that the 0% group has a shear strength between 1.628 and 4.688 higher than the 10% group.
Step 5: Assumptions
Normality
We would like to evaluate the distribution of each treatment.
#NOTE: this function is qqPlot from the car package
qqPlot(tindata$STRENGTH, groups=tindata$METHOD)
qqPlot(tindata$STRENGTH, groups=tindata$ANTIMONY)
- The QQ plots for each of the treatments within each of the factors do not lie closely enough to the theoretical line. The normality assumption is not reasonably met.
Equal Variance
Because we did not determine the normality assumption was reasonably met, we will use Levene’s test of homogeneity.
#this function is in the stats package
leveneTest(STRENGTH ~ factor(ANTIMONY), data = tindata,center=median)leveneTest(STRENGTH ~ METHOD, data = tindata,center=median)Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.- The null hypothesis is that the variance of each treatment is equal. With a p-value this large, we determine the equal variance assumption is met for both main effects.
Complete Example: Paper Towels
Three brands of paper towels were compared to test their strength. A random sample of 26 towels from each of the three brands were tested. Weights were added to the center of the towel as it was held tight. Then drops of water were added. There were three levels of water (0 drops, 5 drops, and 15 drops). The strength was recorded as how much weight the towel held before breaking.
Step 1: Collect the Data and Determine Experimental Design
Consider the design of the data collection: For the One and Two Factor ANOVA we are covering the analysis on completely randomized designs. We want to ensure our data meet that criteria and define the other components of the experimental design.
toweldata<-read.csv("https://raw.githubusercontent.com/kvaranyak4/STAT3220/main/C4%20PaperTowel2.csv",header=T)
names(toweldata)[1] "Brand" "Water" "Strength"head(toweldata)table(toweldata$Brand,toweldata$Water)
0 5 15
Bounty 26 26 26
Comfort 26 26 26
Decorator 26 26 26Experimental design
This is a completely randomized design because experimental units are randomly chosen from each treatment. There two factors, brand, with three levels: Bounty, Comfort, and Decorator and water with three levels 0, 5, and 15 drops. Therefore, there are 9 treatments. It is balanced because each treatment has the same number of experimental units. The response variable is the strength of the towel.
Step 2: Hypothesize Relationship (Exploratory Data Analysis)
For two factor designs, we will examine the interaction plot first. We can look at side by side box plots of the main effect treatment groups and compare the mean value.
# We want to begin with examining the interaction
interaction.plot(toweldata$Brand, toweldata$Water, toweldata$Strength,fun=mean,trace.label="Water Drops"
, xlab="Brand",ylab="Mean Strength",main="Interaction Plot Water Drops X Brand")
interaction.plot(toweldata$Water, toweldata$Brand, toweldata$Strength,fun=mean,trace.label="Brand"
, xlab="Water",ylab="Mean Strength",main="Interaction Plot Water Drops X Brand")
# We can also look at the main effects
#cooling method
boxplot(Strength~Brand, data=toweldata)
tapply(toweldata$Strength,toweldata$Brand,summary)$Bounty
Min. 1st Qu. Median Mean 3rd Qu. Max.
1450 2050 2700 2599 3088 3600
$Comfort
Min. 1st Qu. Median Mean 3rd Qu. Max.
300.0 456.2 1700.0 1772.8 2900.0 3700.0
$Decorator
Min. 1st Qu. Median Mean 3rd Qu. Max.
300 500 700 1123 2000 2700 #antimony level
boxplot(Strength~Water, data=toweldata)
tapply(toweldata$Strength,toweldata$Water,summary)$`0`
Min. 1st Qu. Median Mean 3rd Qu. Max.
1700 2425 2900 2824 3300 3700
$`5`
Min. 1st Qu. Median Mean 3rd Qu. Max.
550 800 1700 1764 2650 3450
$`15`
Min. 1st Qu. Median Mean 3rd Qu. Max.
300.0 400.0 487.5 906.1 1687.5 2350.0 Step 3: Perform the ANOVA
Significance of Interaction
NOTE: With a two factor factorial design, we will always test the significance of the interaction first.
- Null Hypothesis \(H_0 : \mu_{B,0} = \mu_{D,0} = \mu_{C,0} =...=\mu_{C,15}\) (there are 9 means!)
- Alternative Hypothesis at least two means are different
# We include the interaction
# Because water was read in as a numeric variable, we want to
# indicate that it is actually a factor
towelanova1<-aov(Strength~factor(Water)*Brand, data=toweldata)
summary(towelanova1) Df Sum Sq Mean Sq F value Pr(>F)
factor(Water) 2 144038884 72019442 1277.0 <2e-16 ***
Brand 2 85291704 42645852 756.2 <2e-16 ***
factor(Water):Brand 4 27105924 6776481 120.2 <2e-16 ***
Residuals 225 12689471 56398
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Distribution of Test Statistic: F with (3-1)(3-1), 3*3(26-1)= 4,225 DF
- Test Statistic F test: F = 120.2 with pvalue <0.00001
- CONCLUSION: With a pvalue this small, we reject the null hypothesis. There is evidence of an interaction. This means that the strength for each level of brand does depend on the number of drops.
- It is appropriate to perform the post hoc analysis because the design is balanced and the interaction was significant.
Step 4: Post Hoc Analysis
TukeyHSD(towelanova1,conf.level = .95) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Strength ~ factor(Water) * Brand, data = toweldata)
$`factor(Water)`
diff lwr upr p adj
5-0 -1059.9359 -1149.6554 -970.2164 0
15-0 -1918.2692 -2007.9887 -1828.5497 0
15-5 -858.3333 -948.0528 -768.6138 0
$Brand
diff lwr upr p adj
Comfort-Bounty -825.9615 -915.6810 -736.2420 0
Decorator-Bounty -1475.3205 -1565.0400 -1385.6010 0
Decorator-Comfort -649.3590 -739.0785 -559.6395 0
$`factor(Water):Brand`
diff lwr upr p adj
5:Bounty-0:Bounty -171.15385 -377.46198 35.15428 0.1930999
15:Bounty-0:Bounty -1176.92308 -1383.23121 -970.61495 0.0000000
0:Comfort-0:Bounty 157.69231 -48.61582 364.00044 0.2922960
5:Comfort-0:Bounty -1336.53846 -1542.84659 -1130.23033 0.0000000
15:Comfort-0:Bounty -2647.11538 -2853.42351 -2440.80725 0.0000000
0:Decorator-0:Bounty -828.84615 -1035.15428 -622.53802 0.0000000
5:Decorator-0:Bounty -2343.26923 -2549.57736 -2136.96110 0.0000000
15:Decorator-0:Bounty -2601.92308 -2808.23121 -2395.61495 0.0000000
15:Bounty-5:Bounty -1005.76923 -1212.07736 -799.46110 0.0000000
0:Comfort-5:Bounty 328.84615 122.53802 535.15428 0.0000416
5:Comfort-5:Bounty -1165.38462 -1371.69275 -959.07649 0.0000000
15:Comfort-5:Bounty -2475.96154 -2682.26967 -2269.65341 0.0000000
0:Decorator-5:Bounty -657.69231 -864.00044 -451.38418 0.0000000
5:Decorator-5:Bounty -2172.11538 -2378.42351 -1965.80725 0.0000000
15:Decorator-5:Bounty -2430.76923 -2637.07736 -2224.46110 0.0000000
0:Comfort-15:Bounty 1334.61538 1128.30725 1540.92351 0.0000000
5:Comfort-15:Bounty -159.61538 -365.92351 46.69275 0.2765142
15:Comfort-15:Bounty -1470.19231 -1676.50044 -1263.88418 0.0000000
0:Decorator-15:Bounty 348.07692 141.76879 554.38505 0.0000105
5:Decorator-15:Bounty -1166.34615 -1372.65428 -960.03802 0.0000000
15:Decorator-15:Bounty -1425.00000 -1631.30813 -1218.69187 0.0000000
5:Comfort-0:Comfort -1494.23077 -1700.53890 -1287.92264 0.0000000
15:Comfort-0:Comfort -2804.80769 -3011.11582 -2598.49956 0.0000000
0:Decorator-0:Comfort -986.53846 -1192.84659 -780.23033 0.0000000
5:Decorator-0:Comfort -2500.96154 -2707.26967 -2294.65341 0.0000000
15:Decorator-0:Comfort -2759.61538 -2965.92351 -2553.30725 0.0000000
15:Comfort-5:Comfort -1310.57692 -1516.88505 -1104.26879 0.0000000
0:Decorator-5:Comfort 507.69231 301.38418 714.00044 0.0000000
5:Decorator-5:Comfort -1006.73077 -1213.03890 -800.42264 0.0000000
15:Decorator-5:Comfort -1265.38462 -1471.69275 -1059.07649 0.0000000
0:Decorator-15:Comfort 1818.26923 1611.96110 2024.57736 0.0000000
5:Decorator-15:Comfort 303.84615 97.53802 510.15428 0.0002268
15:Decorator-15:Comfort 45.19231 -161.11582 251.50044 0.9989206
5:Decorator-0:Decorator -1514.42308 -1720.73121 -1308.11495 0.0000000
15:Decorator-0:Decorator -1773.07692 -1979.38505 -1566.76879 0.0000000
15:Decorator-5:Decorator -258.65385 -464.96198 -52.34572 0.0035990plot(TukeyHSD(towelanova1,conf.level = .95))
# Compute Minimum Significant Difference w for interactions
# q(p,v) Critical studentized range
q<-qtukey(.95,nmeans=9,df=225)
q[1] 4.429682s<-sqrt(56398)
n<-26
w<-q*(s/sqrt(n))
w[1] 206.3088- We can see that most pairs are significantly different. Whether we use w, condfidence intervals, or p-values. Each approach will lead to the same conclusion at the same significance level.
Step 5: Assumptions
Normality
We would like to evaluate the distribution of each treatment.
# NOTE: this function is qqPlot from the car package
qqPlot(toweldata$Strength, groups=toweldata$Brand)
qqPlot(toweldata$Strength, groups=toweldata$Water)
## EXTRA PLOTTING NOTES
# There are a few options to create plots for each treatment
# So far we have been using base R plotting.
# We can use a loop to make histograms for each level of one factor
par(mfcol=c(2, 2))
for (i in unique(toweldata$Brand)){
hist(toweldata$Strength[toweldata$Brand==i], xlab="Strength", main=paste("Histogram of Brand", i))
}
unique(toweldata$Brand)[1] "Decorator" "Comfort" "Bounty" 
# We can use dev.off() to clear the plot formatting
# dev.off()
# Making plots by two factors is even more complicated in base R plotting
# The package ggplot2 with the function ggplot() is useful for
# more customizable graphics.
#We can do this for histograms
towelhist<-ggplot(toweldata, aes(x=Strength))+ geom_histogram(bins=15, center=.5)+ facet_grid(Brand~Water)
towelhist
# We can do this for QQ plots
towelQQ<-ggplot(toweldata, aes(sample=Strength))+ geom_qq() +
geom_qq_line(color="red")+ facet_grid(Brand~Water)
towelQQ
- The population distributions of the response do not appear to meet the normality assumption.
Equal Variance
Since Normality does not seem reasonably met, we will use Levene’s Test of equal variance.
leveneTest(Strength ~ Brand*factor(Water), data = toweldata,center=median)boxplot(Strength ~ Brand+Water, data = toweldata)
- There is not evidence that the constant variance assumption is met. While these data meet the threshold of treatment sizes of “20”, it is recommended that we either transform the response variable (log, square root, etc) to better fit the assumptions. However, recall when we interpret the response, we cannot give as practical interpretations of the mean of the response. So, instead we may pursue a different “non-parametric” type of analysis.