STAT 3220 Unit 2
These data were collected from 2007 - 2009 by Dr. Kristen Gorman with the Palmer Station Long Term Ecological Research Program, part of the US Long Term Ecological Research Network. The data were imported directly from the Environmental Data Initiative (EDI) Data Portal, and are available for use by CC0 license (“No Rights Reserved”) in accordance with the Palmer Station Data Policy.
Step 1: Collect the Data
Check the appropriateness of response variable for regression: View a histogram of response variable. It should be continuous, and approximately unimodal and symmetric, with few outliers.
pendata<-read.csv("https://raw.githubusercontent.com/kvaranyak4/STAT3220/main/penguins.csv")
head(pendata)names(pendata)## [1] "X" "species" "island"
## [4] "bill_length_mm" "bill_depth_mm" "flipper_length_mm"
## [7] "body_mass_g" "sex" "year"hist(pendata$body_mass_g, xlab="Body Mass", main="Histogram of Body Mass (in grams)") 
Step 2: Hypothesize Relationship (Exploratory Data Analysis)
We explore the box plots and means for each qualitative variable explanatory variable then classify the relationships as existent or not.
#Summary Statistics for response variable grouped by each level of the response
tapply(pendata$body_mass_g,pendata$species,summary)$Adelie
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
2850 3350 3700 3701 4000 4775 1
$Chinstrap
Min. 1st Qu. Median Mean 3rd Qu. Max.
2700 3488 3700 3733 3950 4800
$Gentoo
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
3950 4700 5000 5076 5500 6300 1 tapply(pendata$body_mass_g,pendata$sex,summary)$female
Min. 1st Qu. Median Mean 3rd Qu. Max.
2700 3350 3650 3862 4550 5200
$male
Min. 1st Qu. Median Mean 3rd Qu. Max.
3250 3900 4300 4546 5312 6300 #Box plots for species and sex
boxplot(body_mass_g~species,pendata, ylab="Body Mass (grams)")
boxplot(body_mass_g~sex,pendata, ylab="Body Mass (grams)")
#Scatter plots for quantitative variables
for (i in names(pendata)[4:6]) {
plot(pendata[,i], pendata$body_mass_g,xlab=i,ylab="Body Mass (grams)")
}
#Correlations for quantitative variables
round(cor(pendata[4:6],pendata$body_mass_g,use="complete.obs"),3) [,1]
bill_length_mm 0.595
bill_depth_mm -0.472
flipper_length_mm 0.871There appears to be relationships with each of the qualitative explanatory variables because the mean value of y is different for each level.
bill length and flipper length have positive linear relationships with the response.
bill depth has a negative relationship. however it appears there might be groups in within with positive relationship (we will look at this later)
For the sake of this example, we will hold off on adding the quantitative variables right now.
Therefore, our hypothesized model is:
\(E(bodymass)=\beta_0+\beta_1 SpeciesC+\beta_2 SpeciesG+\beta_3 Sex\)
- where SpeciesC= 1 if Chinstrap, 0 otherwise
- SpeciesG - 1 if Gentoo, 0 otherwise
- Sex = 1 if male, 0 otherwise
Step 3: Estimate the model parameters (fit the model using R)
penmod1<-lm(body_mass_g~species+sex,data=pendata)
summary(penmod1)
Call:
lm(formula = body_mass_g ~ species + sex, data = pendata)
Residuals:
Min 1Q Median 3Q Max
-816.87 -217.80 -16.87 227.61 882.20
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3372.39 31.43 107.308 <2e-16 ***
speciesChinstrap 26.92 46.48 0.579 0.563
speciesGentoo 1377.86 39.10 35.236 <2e-16 ***
sexmale 667.56 34.70 19.236 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 316.6 on 329 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8468, Adjusted R-squared: 0.8454
F-statistic: 606.1 on 3 and 329 DF, p-value: < 2.2e-16The prediction equation is:\(\widehat{bodymass}=3372.39+26.92SpeciesC+1377.86SpeciesG+667.56Sex\)
- where SpeciesC= 1 if Chinstrap, 0 otherwise
- SpeciesG = 1 if Gentoo, 0 otherwise
- Sex = 1 if male, 0 otherwise
Interpret the estimates
- The average body mass for Adelie male penguins is 3372.39 grams.
- The estimated mean body mass for Chinstrap penguins is 26.92 grams higher than Adelie penguins, given gender is held constant.
- The estimated mean body mass for Gentoo penguins is 1377.86 grams higher than Adelie penguins, given gender is held constant.
- The estimated mean body mass for male penguins is 667.56 grams higher than female penguins, given species is held constant.
Step 4: Specify the distribution of the errors and find the estimate of the variance
- \(\epsilon \overset{\mathrm{iid}}{\sim} N(0,\sigma^2 )\)
- estimate of \(\sigma\) is \(\sqrt{MSE}=316.6\)
- Approximately 95% of our predictions will be within 2*316.6=633.2 grams of the actual body mass for the penguins.
Step 5: Evaluate the Utility of the model
First we Perform the Global F Test:
Hypotheses:
- \(H_0: \beta_1= \beta_2=\beta_3=0\) (the model is not adequate)
- \(H_a\):at least one of \(\beta_1 , \beta_2 , \beta_3\) does not equal 0 (the model is adequate)
Distribution of test statistic: F with 3,329 DF
Test Statistic: F=606.1
Pvalue: <2.2e-16
Decision: pvalue<0.05 -> REJECT H0
Conclusion: The model with species and sex is adequate at body mass for the penguins.
What terms can we test “indivudally?”
Species is determined by two parameters (SpeciesC and SpeciesG dummy variables), so we would want to keep or remove BOTH of those dummy variables. We will cover this later.
We can test sex because it is only one parameter in the model.
Hypotheses:
- \(H_0: \beta_3=0\)
- \(H_a:\beta_3 \neq 0\)
Distribution of test statistic: T with 329 DF
Test Statistic: t=19.236
Pvalue: <2.2e-16
Decision: pvalue<0.05 -> REJECT H0
Conclusion: The sex of the penguin is significant at predicting the body mass of a penguin given species is a constant in the model. We will keep it in the model.
Quantative predictors
Note we can include both qualitative and quantitative predictors in the model.
penmod2<-lm(body_mass_g~species+sex+bill_depth_mm,data=pendata)
summary(penmod2)
Call:
lm(formula = body_mass_g ~ species + sex + bill_depth_mm, data = pendata)
Residuals:
Min 1Q Median 3Q Max
-851.2 -187.7 -3.5 221.7 936.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1576.99 357.08 4.416 1.37e-05 ***
speciesChinstrap 19.44 44.87 0.433 0.665
speciesGentoo 1721.70 77.88 22.106 < 2e-16 ***
sexmale 513.99 45.24 11.360 < 2e-16 ***
bill_depth_mm 102.04 20.22 5.046 7.48e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 305.5 on 328 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8578, Adjusted R-squared: 0.8561
F-statistic: 494.8 on 4 and 328 DF, p-value: < 2.2e-16Interpretation: For each millimeter increase in bill depth, we expect the body weight of a penguin to increase by 102.04 grams, given species and sex are held constant.
plot(body_mass_g~bill_depth_mm , col=factor(species),data=pendata,xlab="Bill Depth",ylab="Body Mass (grams)")
legend("topright",legend = levels(factor(pendata$species)), pch = 19,
col = factor(levels(factor(pendata$species))))
NOTE: When we regroup the scatter plot for Bill Depth by species, the relationships for each level are positive. Note also we see each species would have a similar slope, but different y-intercept.
Let’s suppose we simplify the model to just species and bill depth:
\(E(body mass)=\beta_0+\beta_1SpeciesC+ \beta_2SpeciesG+\beta_3 BillDepth\) , where SpeciesC= 1 if Chinstrap, 0 otherwise and SpeciesG = 1 if Gentoo, 0 otherwise
- The body mass- bill depth y intercept for Adelie (SpeciesC=0,SpeciesG=0) will be: \(\beta_0\)
- The body mass- bill depth y intercept for Chinstrap (SpeciesC=1,SpeciesG=0) will be: \(\beta_0+\beta_1\)
- The body mass- bill depth y intercept for Gentroo (SpeciesC=0,SpeciesG=1) will be: \(\beta_0+\beta_2\)
- All species will have the same body mass- bill depth slope of \(beta_3\)
What if there were different slopes for each group? We will explore this next time!