STAT 3220 Unit 3
Dr. Kristen Gorman with the Palmer Station Long Term Ecological Research Program set out to study the impact of Ecological Sexual Dimorphism in arctic penguins. In 2014, that particular study compared deep-dive foraging habits of male and female penguins and how that resulted in different biological outcomes. To read more abut it, you can go to read the full study here. And you can find more of Dr. Gorman’s more recent work here
We only have a portion of the data from the larger scale study, so our goals and scope of analysis will be different. Our aim is to determine which factors impact the body mass of a penguin. This is useful to know as body mass is an important indicator of penguin’s health.
Our data table contains the following columns:
- Observation: the observation number in the data table
- body_mass_g: an integer denoting the weight of a penguin’s body (grams) (Response Variable)
- species: a factor denoting the penguin species (Adélie, Chinstrap, or Gentoo)
- island: a factor denoting the Palmer Archipelago island in Antarctica where each penguin was observed (Biscoe Point, Dream Island, or Torgersen Island)
- bill_length_mm: a number denoting length of the dorsal ridge of a penguin bill (millimeters)
- bill_depth_mm: a number denoting the depth of a penguin bill (millimeters)
- flipper_length_mm: an integer denoting the length of a penguin flipper (millimeters)
- sex: a factor denoting the sex of a penguin sex (male, female) based on molecular data
- year: an integer denoting the year of study (2007, 2008, or 2009)
Step 1-5 Recap
pendata<-read.csv("https://raw.githubusercontent.com/kvaranyak4/STAT3220/main/penguins.csv")
head(pendata)names(pendata)[1] "Observation" "body_mass_g" "species"
[4] "island" "bill_length_mm" "bill_depth_mm"
[7] "flipper_length_mm" "sex" "year" Hypothesized Linear Model After going through our model building process:
\(E(bodymass)=\beta_0+\beta_1BillDep+\beta_2BillLen+\beta_3FlipperLen+\beta_4SpeciesC+\beta_5SpeciesG+\beta_6SexM+\beta_7BillDep*BillLen\)
- where Species C= 1 if Chinstrap, 0 otherwise
- SpeciesG = 1 if Gentoo, 0 otherwise
- SexM = 1 if male, 0 otherwise
Step 6: Check Model Assumptions (Residual Analysis)
Unit 3.1
Check the model assumptions.
#Fit the final model
penmodfinal<-lm(body_mass_g~bill_depth_mm+bill_length_mm+flipper_length_mm+bill_depth_mm*bill_length_mm+species+sex,data=pendata)
#There are a few options to view plots for assumptions
#Residuals Plots of explanatory variables vs residuals
residualPlots(penmodfinal,tests=F)
#Residual vs Fitted and QQ plot
plot(penmodfinal, which=c(1,2))
#histogram of residuals
hist(residuals(penmodfinal))
Lack of Fit - Plots to evaluate: Residual Plots - There were no clear trends that could be modeled in the residual plots. The mean 0 (lack of fit) assumption was not violated.
Constant Variance - Plots to evaluate: Residual Plots and Residual vs Fitted - There were no clear fanning patterns in the residual plots nor the residual vs predicted plots. For qualitative variables we look for even spread of residuals across the levels. The constant variance (homoscedasticity) assumption was not violated.
Normality - Plots to evaluate: QQ Plot and Histogram of Residuals - The observations fall linearly on the QQ plot. (The residuals are unimodal and symmetric in the histogram). The normality assumption is not violated.
Independence - Plots to evaluate: No specific plots. We could perform the Durbin Watson test, but we will not do that in the scope of this course. - We do not believe the information of one penguin depends on another penguin. The independence assumption is likely not violated.
Notes on Corrections Violations of Assumptions
In this example there do not appear to be severe violations of the assumptions. Therefore, we do not need to transform the response or explanaotry variables We can proceed to assessing extreme observations.
Unit 3.2
Identify outliers and influential observations.
Label the observations that exceed the threshold for each statistic:
- Cook’s Distance: Identifies influential
observations. Threshold is 4/n
- Circles on “influence plot” are indicative of the Cook’s distance value
- Leverage (Hat): Identifies outliers in the x direction (extreme values of explanatory variables). Threshold is 2(k+1)/n
- Studentizied Residuals: Identifies outliers in the y direction. Threshold is +/- 2
- Deleted Studentizied Residuals (R Student): Identifies general outliers and influential observations. Threshold is +/- 2
Computing the Thresholds
n will be how many observations were used to fit the model, not necessarily how many observations were in the original data table. To determine this number we can go to the model F degrees of freedom k, n-k-1. In the pen model we have 7, 325. So k=7 and n=333.
#Determine n
summary(penmodfinal)
Call:
lm(formula = body_mass_g ~ bill_depth_mm + bill_length_mm + flipper_length_mm +
bill_depth_mm * bill_length_mm + species + sex, data = pendata)
Residuals:
Min 1Q Median 3Q Max
-706.90 -175.12 -6.14 176.80 880.08
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6290.446 1832.633 -3.432 0.000675 ***
bill_depth_mm 333.924 98.231 3.399 0.000760 ***
bill_length_mm 130.274 41.058 3.173 0.001653 **
flipper_length_mm 16.535 2.888 5.725 2.35e-08 ***
speciesChinstrap -200.763 82.322 -2.439 0.015273 *
speciesGentoo 923.788 132.379 6.978 1.68e-11 ***
sexmale 391.603 47.370 8.267 3.55e-15 ***
bill_depth_mm:bill_length_mm -6.344 2.290 -2.770 0.005919 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 284.4 on 325 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8778, Adjusted R-squared: 0.8752
F-statistic: 333.7 on 7 and 325 DF, p-value: < 2.2e-16#numDF + (k+1) = n-k-1+(k+1) = n
n<-df.residual(penmodfinal)+length(coef(penmodfinal))
n[1] 333#Leverage Cut off
levcut<-(2*(length(coef(penmodfinal))))/n
levcut[1] 0.04804805#Cook's Distance Cutoff
cdcut<-4/n
cdcut[1] 0.01201201Now we can locate these values on the influence plots to identify outliers and influential observations.
#Cooks Distance Thresholds
plot(penmodfinal,which=4)
#Leverage vs Studentized Residuals
influencePlot(penmodfinal,fill=F)
#Deleted Studentized Residals
ols_plot_resid_stud_fit(penmodfinal)
Based on the labeled observations from the plots:
With respect to Leverage, observations 314, 15, 295 (and others) are considered “leverage”=outliers in x because their values exceed the threshold of 2(7+1)/333=0.048.
With respect to Studentizied Residuals, observations 170, 40, 315, 314 (and others) are considered outliers in y because their values exceed the threshold of +/-2 .
With respect to Cook’s Distance, observations 293, 314, 315 are considered influential because their values exceed the threshold of 4/n=4/333=0.012.
With respect to Deleted Studentizied Residuals, observations that are labelled red (164, 35, 304, etc) are considered influential because their values exceed the threshold of +/-2.
Note: The “influence output” only labels some of the observations that exceed our thresholds.
Determine all of the observations that exceed our thresholds
We can use the which() function to print values that
meet certain criteria from a logical expression.
#General expressions about the variables in the data table
#Quantitative Variables
which(pendata$bill_depth_mm<17.3)
# Qualitative Variables
which(pendata$species=="Gentoo")
#Multiple arguments
which(pendata$species=="Gentoo"&pendata$bill_depth_mm<17.3)
which(pendata$species=="Gentoo"|pendata$bill_depth_mm<17.3)We have several functions that will compute the values of our influence diagnostics for every observation. The input argument for each of these functions is the model that was fit.
- Leverage:
hatvalues() - Studentized Residuals:
stdres() - Cook’s Distance:
cooks.distance() - Deleted Studentized Residuals:
rstudent()
# Using our influential diagnostics
## Leverage
which(hatvalues(penmodfinal)>16/333) 15 20 115 130 143 186 190 240 288 294 308 314 340
10 15 109 124 137 179 183 232 277 283 297 303 329 ## Studentized Residual
which(abs(stdres(penmodfinal))>2) 8 40 46 82 86 94 105 110 120 122 132 166 170 196 232 285 293 314 315 325
7 35 41 76 80 88 99 104 114 116 126 160 164 189 224 274 282 303 304 314 ## Cooks Distance
which(cooks.distance(penmodfinal)>4/333) 40 82 112 118 122 166 170 193 232 285 293 294 297 298 314 315 325
35 76 106 112 116 160 164 186 224 274 282 283 286 287 303 304 314 ## Deleted Studentized Residual
which(abs(rstudent(penmodfinal))>2) 8 40 46 82 86 94 105 110 120 122 132 166 170 196 232 285 293 314 315 325
7 35 41 76 80 88 99 104 114 116 126 160 164 189 224 274 282 303 304 314 #Multiple arguments
which(hatvalues(penmodfinal)>16/333&cooks.distance(penmodfinal)>4/333&abs(rstudent(penmodfinal))>2&abs(stdres(penmodfinal))>2)314
303 The top row is the original row position of the observation from the data frame, the bottom row is the observation’s position index on the list of influence diagnostics. These do not match in this example because there are NA values in the data frame.
Observations 15, 20, 115, 130, 143, 186, 190, 240, 288, 294, 308, 314, 340 exceed the threshold for leverage.
Observations 40, 82, 112, 118, 122, 166, 170, 193, 232, 285, 293, 294, 297, 298, 314, 315, 325 exceed the threshold for Cook’s Distance.
Observations 8, 40, 46, 82, 86, 94, 105, 110, 120, 122, 132, 166, 170, 196, 232, 285, 293, 314, 315, 325 exceed the threshold for Deleted Studentized Residuals.
Observations 8, 40, 46, 82, 86, 94, 105, 110, 120, 122, 132, 166, 170, 196, 232, 285, 293, 314, 315, 325 exceed the threshold for Studentized Residuals.
Observation 314 exceeds in all 4 categories.
Evaluate Extreme Observations Closer
Next steps would be to evaluate closer with dffits or dfbetas and see how those observations compare to the other observations in the data. Then look at the EDA. Determine if the observations are very influential and/or truly representative of the data. We should NOT just remove without further exploration.
#First we will store the DFFits to make them easier to work with
pendffits<-dffits(penmodfinal)
# Calling the observation numbers (row names)
# don't get too worried on the syntax
pendffits[c("15","40","170","294","314","315")] 15 40 170 294 314 315
0.1286874 0.3818284 0.3833745 0.3424947 0.5436063 -0.4620501 # Now we can compare these values to the rest
summary(dffits(penmodfinal)) Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.462050 -0.093274 -0.003084 0.001633 0.101550 0.543606 # Slightly different syntax for DFBetas
# Because DFBetas is a matrix. not just a vector
pendfbeta<-dfbetas(penmodfinal)
pendfbeta[c("15","40","170","294","314","315"),] (Intercept) bill_depth_mm bill_length_mm flipper_length_mm speciesChinstrap
15 -0.08557838 0.09310769 0.06910227 0.02660027 0.0655398436
40 0.02653293 0.05209415 0.03999337 -0.25112153 -0.0043988499
170 -0.12168270 0.11147375 0.13543939 0.04868884 0.0246383890
294 -0.01852371 0.02920232 0.08105004 -0.17108150 -0.1584788882
314 0.26622725 -0.32579551 -0.36372431 0.14860592 -0.0004064124
315 0.21521173 -0.21042438 -0.25248533 -0.04511119 -0.1638646983
speciesGentoo sexmale bill_depth_mm:bill_length_mm
15 0.02492217 0.009146321 -0.08271890
40 0.12938768 0.170995071 -0.04103063
170 -0.05802485 0.169876311 -0.13743246
294 -0.04443237 -0.119146477 -0.03129413
314 0.07753382 -0.100034429 0.36711794
315 0.18047644 -0.001915482 0.25222563summary(dfbetas(penmodfinal)) (Intercept) bill_depth_mm bill_length_mm
Min. :-2.241e-01 Min. :-0.3257955 Min. :-3.637e-01
1st Qu.:-1.834e-02 1st Qu.:-0.0121734 1st Qu.:-9.957e-03
Median :-3.307e-04 Median : 0.0009279 Median : 2.128e-04
Mean : 7.712e-05 Mean :-0.0000870 Mean :-7.541e-05
3rd Qu.: 1.552e-02 3rd Qu.: 0.0194839 3rd Qu.: 1.921e-02
Max. : 2.662e-01 Max. : 0.2017456 Max. : 2.407e-01
flipper_length_mm speciesChinstrap speciesGentoo
Min. :-2.511e-01 Min. :-2.327e-01 Min. :-3.006e-01
1st Qu.:-1.948e-02 1st Qu.:-1.623e-02 1st Qu.:-1.962e-02
Median : 1.204e-04 Median :-7.051e-05 Median : 6.127e-04
Mean :-3.359e-05 Mean :-9.018e-05 Mean : 8.558e-06
3rd Qu.: 2.059e-02 3rd Qu.: 1.830e-02 3rd Qu.: 1.995e-02
Max. : 2.652e-01 Max. : 2.515e-01 Max. : 2.159e-01
sexmale bill_depth_mm:bill_length_mm
Min. :-2.001e-01 Min. :-2.387e-01
1st Qu.:-2.291e-02 1st Qu.:-1.882e-02
Median :-1.235e-03 Median :-4.818e-04
Mean :-9.972e-05 Mean : 9.593e-05
3rd Qu.: 1.990e-02 3rd Qu.: 1.117e-02
Max. : 2.072e-01 Max. : 3.671e-01 # Look at those values in the table
pendata[c("15","40","170","294","314","315"),]summary(pendata) Observation body_mass_g species island
Min. : 1.00 Min. :2700 Length :344 Length :344
1st Qu.: 86.75 1st Qu.:3550 N.unique : 3 N.unique : 3
Median :172.50 Median :4050 N.blank : 0 N.blank : 0
Mean :172.50 Mean :4202 Min.nchar: 6 Min.nchar: 5
3rd Qu.:258.25 3rd Qu.:4750 Max.nchar: 9 Max.nchar: 9
Max. :344.00 Max. :6300
NAs :2
bill_length_mm bill_depth_mm flipper_length_mm sex
Min. :32.10 Min. :13.10 Min. :172.0 Length :344
1st Qu.:39.23 1st Qu.:15.60 1st Qu.:190.0 N.unique : 2
Median :44.45 Median :17.30 Median :197.0 N.blank : 0
Mean :43.92 Mean :17.15 Mean :200.9 Min.nchar: 4
3rd Qu.:48.50 3rd Qu.:18.70 3rd Qu.:213.0 Max.nchar: 6
Max. :59.60 Max. :21.50 Max. :231.0 NAs : 11
NAs :2 NAs :2 NAs :2
year
Min. :2007
1st Qu.:2007
Median :2008
Mean :2008
3rd Qu.:2009
Max. :2009
Now we can compare the DFFits and DFBetas to the rest. Additionally, we can compare the x values to the rest of the data. Upon investigating, observation 314 appears to be very influential and has high values of DFFits. We may consider removing this.
Subsetting our data frame to remove influential observations
We can use the which() function inside our square
bracket notation to create new data frames based on subsetting
requirements. We can store these as new data objects, or use the
subsetting directly into another function’s arguments.
# General formatting for subsetting with the which() function
## Quantitative Variables
quantsubsetpen<-pendata[which(pendata$bill_depth_mm<17.3),]
#Example
plot((body_mass_g)~flipper_length_mm,data=quantsubsetpen)
plot((body_mass_g)~flipper_length_mm,data=pendata[which(pendata$bill_depth_mm<17.3),])
# Qualitative Variables
qualsubsetpen<-pendata[which(pendata$species=="Gentoo"),]
# Multiple Arguments
multaubsetpen1<-pendata[which(pendata$species=="Gentoo"&pendata$bill_depth_mm<17.3),]
multaubsetpen2<-pendata[which(pendata$species=="Gentoo"|pendata$bill_depth_mm<17.3),]If you want to remove just one or more observations, or a set of
observations, we can use a - (“minus”) or the which()
function within a subset.
# You can do this directly into the data argument in lm
# remove one single value, based on the observation's original row number
penmodrem.1<-lm(body_mass_g~bill_depth_mm+bill_length_mm+flipper_length_mm+bill_depth_mm*bill_length_mm+species+sex,data=pendata[-314,])
summary(penmodrem.1)
Call:
lm(formula = body_mass_g ~ bill_depth_mm + bill_length_mm + flipper_length_mm +
bill_depth_mm * bill_length_mm + species + sex, data = pendata[-314,
])
Residuals:
Min 1Q Median 3Q Max
-704.2 -170.2 0.9 175.3 875.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6775.864 1838.196 -3.686 0.000267 ***
bill_depth_mm 365.764 98.923 3.697 0.000256 ***
bill_length_mm 145.132 41.470 3.500 0.000531 ***
flipper_length_mm 16.108 2.881 5.592 4.78e-08 ***
speciesChinstrap -200.730 81.904 -2.451 0.014782 *
speciesGentoo 913.577 131.798 6.932 2.26e-11 ***
sexmale 396.318 47.184 8.399 1.42e-15 ***
bill_depth_mm:bill_length_mm -7.180 2.313 -3.104 0.002079 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 283 on 324 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8793, Adjusted R-squared: 0.8767
F-statistic: 337.1 on 7 and 324 DF, p-value: < 2.2e-16# note DF is different
#remove more than one, based on the observations' original row numbers
penmodrem.2<-lm(body_mass_g~bill_depth_mm+bill_length_mm+flipper_length_mm+bill_depth_mm*bill_length_mm+species+sex,data=pendata[-c(15,40,170,294,314,315),])
summary(penmodrem.2)
Call:
lm(formula = body_mass_g ~ bill_depth_mm + bill_length_mm + flipper_length_mm +
bill_depth_mm * bill_length_mm + species + sex, data = pendata[-c(15,
40, 170, 294, 314, 315), ])
Residuals:
Min 1Q Median 3Q Max
-681.17 -171.80 1.54 171.99 688.00
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6820.422 1833.965 -3.719 0.000236 ***
bill_depth_mm 359.304 99.122 3.625 0.000336 ***
bill_length_mm 142.566 41.217 3.459 0.000616 ***
flipper_length_mm 17.212 2.859 6.020 4.79e-09 ***
speciesChinstrap -181.755 81.946 -2.218 0.027260 *
speciesGentoo 883.023 128.661 6.863 3.52e-11 ***
sexmale 385.230 46.275 8.325 2.51e-15 ***
bill_depth_mm:bill_length_mm -7.107 2.309 -3.079 0.002260 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 274.2 on 319 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8845, Adjusted R-squared: 0.882
F-statistic: 349 on 7 and 319 DF, p-value: < 2.2e-16#remove based on logical expression
penmodrem.3<-lm(body_mass_g~bill_depth_mm+bill_length_mm+flipper_length_mm+bill_depth_mm*bill_length_mm+species+sex,data=pendata[which(pendata$bill_depth_mm<17.3),])
summary(penmodrem.3)
Call:
lm(formula = body_mass_g ~ bill_depth_mm + bill_length_mm + flipper_length_mm +
bill_depth_mm * bill_length_mm + species + sex, data = pendata[which(pendata$bill_depth_mm <
17.3), ])
Residuals:
Min 1Q Median 3Q Max
-683.38 -162.91 -14.99 155.11 860.97
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4913.958 4618.410 -1.064 0.28901
bill_depth_mm 254.936 281.470 0.906 0.36650
bill_length_mm 82.344 94.554 0.871 0.38519
flipper_length_mm 17.357 4.703 3.691 0.00031 ***
speciesChinstrap -117.223 135.779 -0.863 0.38930
speciesGentoo 1036.777 213.654 4.853 2.98e-06 ***
sexmale 447.063 71.778 6.228 4.34e-09 ***
bill_depth_mm:bill_length_mm -3.869 6.021 -0.642 0.52152
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 269 on 153 degrees of freedom
(5 observations deleted due to missingness)
Multiple R-squared: 0.917, Adjusted R-squared: 0.9132
F-statistic: 241.4 on 7 and 153 DF, p-value: < 2.2e-16There are TONS of ways we can use R to streamline (or potentially over-complicate) this process. We will stick with removing just a few observations at once.
Cautions on removing observations
- We want to retain as much data as possible
- If we remove some extreme observations then other observations may be identified as extreme
Residual Analysis Conclusions
We may need to iteratively recheck model assumptions and influential diagnostics as we transform variables and/or remove observations. It is good practice to make one change at a time. There is not necessarily one right way to end model building in the residual analysis.
Re-perform the residual analysis without observation 314
#Residuals Plots of explanatory variables vs residuals
residualPlots(penmodrem.1,tests=F)
#Residual vs Fitted and QQ plot
plot(penmodrem.1, which=c(1,2))
#Cooks Distance Thresholds
plot(penmodrem.1,which=4)
#Leverage vs Studentized Residuals
influencePlot(penmodrem.1,fill=F)
#Deleted Studentized Residals
ols_plot_resid_stud_fit(penmodrem.1)
The assumptions don’t seem much different and we still have pretty extreme values. It is likely not worth it to alter this model.
Unit 3.3
Interpretations on Variable Transformations
#Within lm() we can add both transformations to x and or y
#sqrt(x)
penmodTX<-lm(body_mass_g~bill_depth_mm+bill_length_mm+sqrt(flipper_length_mm)+bill_depth_mm*bill_length_mm+species+sex,data=pendata)
summary(penmodTX)
Call:
lm(formula = body_mass_g ~ bill_depth_mm + bill_length_mm + sqrt(flipper_length_mm) +
bill_depth_mm * bill_length_mm + species + sex, data = pendata)
Residuals:
Min 1Q Median 3Q Max
-709.13 -179.20 -5.68 175.75 881.73
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9544.358 2081.033 -4.586 6.45e-06 ***
bill_depth_mm 335.228 98.335 3.409 0.000734 ***
bill_length_mm 131.016 41.104 3.187 0.001575 **
sqrt(flipper_length_mm) 461.331 81.442 5.665 3.25e-08 ***
speciesChinstrap -204.408 82.403 -2.481 0.013622 *
speciesGentoo 931.382 132.169 7.047 1.10e-11 ***
sexmale 392.423 47.405 8.278 3.29e-15 ***
bill_depth_mm:bill_length_mm -6.366 2.292 -2.777 0.005806 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 284.7 on 325 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8776, Adjusted R-squared: 0.875
F-statistic: 332.9 on 7 and 325 DF, p-value: < 2.2e-16#sqrt(y)
penmodTY<-lm(sqrt(body_mass_g)~bill_depth_mm+bill_length_mm+flipper_length_mm+bill_depth_mm*bill_length_mm+species+sex,data=pendata)
summary(penmodTY)
Call:
lm(formula = sqrt(body_mass_g) ~ bill_depth_mm + bill_length_mm +
flipper_length_mm + bill_depth_mm * bill_length_mm + species +
sex, data = pendata)
Residuals:
Min 1Q Median 3Q Max
-6.4098 -1.4023 -0.0304 1.3594 6.0240
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.14762 14.24113 -0.993 0.32124
bill_depth_mm 2.51499 0.76334 3.295 0.00109 **
bill_length_mm 0.95538 0.31906 2.994 0.00296 **
flipper_length_mm 0.12317 0.02244 5.489 8.17e-08 ***
speciesChinstrap -1.41805 0.63971 -2.217 0.02734 *
speciesGentoo 7.28823 1.02870 7.085 8.68e-12 ***
sexmale 3.00334 0.36810 8.159 7.48e-15 ***
bill_depth_mm:bill_length_mm -0.04659 0.01779 -2.618 0.00926 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.21 on 325 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8729, Adjusted R-squared: 0.8701
F-statistic: 318.8 on 7 and 325 DF, p-value: < 2.2e-16#Prediction
predict(penmodTY,interval="predict",level=.95)[1,]Warning in predict.lm(penmodTY, interval = "predict", level = 0.95): predictions on current data refer to _future_ responses fit lwr upr
61.47405 57.07251 65.87559 57.07251^2[1] 3257.27165.87559^2[1] 4339.593#Estimation
predict(penmodTY,interval="confidence",level=.95)[1,] fit lwr upr
61.47405 60.79217 62.15594 - When we transform the explanatory variables, the interpretations do not change drastically
- Coefficient Interpretation: For a one unit increase in the transformed value of X (square root of flipper length), the predicted value of y changes by \(\hat{\beta_i}\) given all else is constant
- For a one unit increase in the square root of flipper length, the predicted value of body mass changes by 461.331 grams, given bill length, bill depth, species and sex are constant.
- Similarly, the transformation for the response variable, y, follow:
- Coefficient Interpretation: For a one unit increase in the value of x, the predicted value of the transformed value of y (y*) changes by \(\hat{\beta_i}\)
- For a one mm increase in the flipper length, the predicted value of the square root of body mass changes by 0.12317 square root of grams, given bill length, bill depth, species and sex are constant.
- Prediction Interval: When we have this set of explanatory variables we expect y* to be between the interval computed.
- We are 95% confident that for an Adelie male penguin with a bill length of 39.1mm, bill depth of 18.7mm, flipper length of 181mm, the square root of its body mass will be between 57.07251 and 65.87559 square root grams.
- Prediction Interval Transformed: When we have this set of explanatory variables we expect the value y to be between the interval computed if we take the inverse transformation of the end points.
- We are 95% confident that for an Adelie male penguin with a bill length of 39.1mm, bill depth of 18.7mm, flipper length of 181mm, its body mass will be between \((57.07251)^2=3257.271\) and \((65.87559)^2=4339.593\) grams.
- Confidence Interval: When we have this set of explanatory variables we expect the average of y* (the transformed value of y) to be between the interval computed.
- We are 95% confident that the true average square root of body mass of Adelie male penguins with a bill length of 39.1mm, bill depth of 18.7mm, flipper length of 181mm, between 60.79217 and 62.15594 square root grams.
- Confidence Interval Transformed: we cannot perform the inverse of the transformation to get the original units.
Notes on Extrapolation
We do not know how the relationship of the variables will behave beyond the experimental region of the data collected.
- Experimental Region: For a single independent variable, the experimental region is the range of the values of x in the sample. For multiple linear regression, the experimental region is the levels of each variable that jointly define the region.
- Extrapolation is when you predict y for values (or estimate for E(y)) using a set of independent variables that are outside the experimental region
- Hidden Extrapolation: x’s fall within the range of each individual x, but outside the jointly defined experimental region
We want to look at the combination of bill length, bill depth, flipper length, species and sex.
#Summarize each variable in the table
summary(pendata) Observation body_mass_g species island
Min. : 1.00 Min. :2700 Length :344 Length :344
1st Qu.: 86.75 1st Qu.:3550 N.unique : 3 N.unique : 3
Median :172.50 Median :4050 N.blank : 0 N.blank : 0
Mean :172.50 Mean :4202 Min.nchar: 6 Min.nchar: 5
3rd Qu.:258.25 3rd Qu.:4750 Max.nchar: 9 Max.nchar: 9
Max. :344.00 Max. :6300
NAs :2
bill_length_mm bill_depth_mm flipper_length_mm sex
Min. :32.10 Min. :13.10 Min. :172.0 Length :344
1st Qu.:39.23 1st Qu.:15.60 1st Qu.:190.0 N.unique : 2
Median :44.45 Median :17.30 Median :197.0 N.blank : 0
Mean :43.92 Mean :17.15 Mean :200.9 Min.nchar: 4
3rd Qu.:48.50 3rd Qu.:18.70 3rd Qu.:213.0 Max.nchar: 6
Max. :59.60 Max. :21.50 Max. :231.0 NAs : 11
NAs :2 NAs :2 NAs :2
year
Min. :2007
1st Qu.:2007
Median :2008
Mean :2008
3rd Qu.:2009
Max. :2009
#Example of Experimental Region for Species and Bill Length
boxplot(bill_length_mm~species,pendata)
tapply(pendata$bill_length_mm,pendata$species,summary)$Adelie
Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
32.10 36.75 38.80 38.79 40.75 46.00 1
$Chinstrap
Min. 1st Qu. Median Mean 3rd Qu. Max.
40.90 46.35 49.55 48.83 51.08 58.00
$Gentoo
Min. 1st Qu. Median Mean 3rd Qu. Max. NAs
40.90 45.30 47.30 47.50 49.55 59.60 1 #Example of Experimental Region for Bill Length and Bill Depth
plot(bill_length_mm~bill_depth_mm,data=pendata)
- We should only use predictions for the range of x values for the data
- Hidden extrapolation example for Species and bill length: We should only predict Adelie penguins with bill lengths between 32.1 mm and 46mm.
- Hidden extrapolation example for bill depth and Bill length: We should only predict penguins that have a bill depth between about 13mm and 15mm with bill lengths from about 40mm to 50mm. But for penguins with a bill depth larger than about 15mm, and range of bill length is reasonable.
Notes on P-values
summary(penmodfinal)
Call:
lm(formula = body_mass_g ~ bill_depth_mm + bill_length_mm + flipper_length_mm +
bill_depth_mm * bill_length_mm + species + sex, data = pendata)
Residuals:
Min 1Q Median 3Q Max
-706.90 -175.12 -6.14 176.80 880.08
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6290.446 1832.633 -3.432 0.000675 ***
bill_depth_mm 333.924 98.231 3.399 0.000760 ***
bill_length_mm 130.274 41.058 3.173 0.001653 **
flipper_length_mm 16.535 2.888 5.725 2.35e-08 ***
speciesChinstrap -200.763 82.322 -2.439 0.015273 *
speciesGentoo 923.788 132.379 6.978 1.68e-11 ***
sexmale 391.603 47.370 8.267 3.55e-15 ***
bill_depth_mm:bill_length_mm -6.344 2.290 -2.770 0.005919 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 284.4 on 325 degrees of freedom
(11 observations deleted due to missingness)
Multiple R-squared: 0.8778, Adjusted R-squared: 0.8752
F-statistic: 333.7 on 7 and 325 DF, p-value: < 2.2e-16The p-value is probability of observing a value of the test statistic at least as contradictory to the null hypothesis as the observed test statistic value, assuming the null hypothesis is true.
- Interpret the p-value for the interaction of bill depth X bill length: The probability of obtaining our test statistic (t=-2.77) or something more extreme (larger in absolute value) from a t distribution with 325 DF if the null is true (\(\beta_7=0\) or the interaction of bill depth and bill length does not contribute to body mass) is 0.0059. This probability is so “small” that it is unlikely the null is true. Therefore we have evidence to support the alternative hypothesis that \(\beta_7 \neq 0\), meaning the interaction does contribute to predicting the response.