── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Residual and Diagnostic Analysis of the Blood Loss Data
We looked at much of this before, but let us go through the above steps systematically. Recall the data set (we want to predict blood loss from weight):
#### Residual and Diagnostic Analysis of the Blood Loss Datadat_thyroid <-read.table(text="weight time blood_loss 44.3 105 503 40.6 80 490 69.0 86 471 43.7 112 505 50.3 109 482 50.2 100 490 35.4 96 513 52.2 120 464", header=TRUE) |>mutate(# create data idsid =1:n() )# show the structure of the data.framestr(dat_thyroid)
'data.frame': 8 obs. of 4 variables:
$ weight : num 44.3 40.6 69 43.7 50.3 50.2 35.4 52.2
$ time : int 105 80 86 112 109 100 96 120
$ blood_loss: int 503 490 471 505 482 490 513 464
$ id : int 1 2 3 4 5 6 7 8
# ggplot: Plot the data with linear regression fit and confidence bandslibrary(ggplot2)p <-ggplot(dat_thyroid, aes(x = weight, y = blood_loss, label = id))p <- p +geom_point()# plot labels next to pointsp <- p +geom_text(hjust =0.5, vjust =-0.5)# plot regression line and confidence bandp <- p +geom_smooth(method = lm)print(p)
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Clearly the heaviest individual is an unusual value that warrants
a closer look (maybe data recording error).
I might be inclined to try a transformation here (such as `log(weight)`) to
make that point a little less influential.
Do any obvious transformations of the data. We will look at transformations later.
Fit the least squares equation. Blood Loss appears significantly negatively associated with weight.
lm_blood_wt <-lm(blood_loss ~ weight, data = dat_thyroid)# use summary() to get t-tests of parameters (slope, intercept)summary(lm_blood_wt)
Call:
lm(formula = blood_loss ~ weight, data = dat_thyroid)
Residuals:
Min 1Q Median 3Q Max
-20.565 -6.189 4.712 8.192 9.382
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 552.4420 21.4409 25.77 2.25e-07 ***
weight -1.3003 0.4364 -2.98 0.0247 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.66 on 6 degrees of freedom
Multiple R-squared: 0.5967, Adjusted R-squared: 0.5295
F-statistic: 8.878 on 1 and 6 DF, p-value: 0.02465
1. *Graphs: Check Standardized Residuals (or the Deleted Residuals).*
The residual plots:
Do you see curvature? There does not appear to be curvature (and it could be hard to detect with so few points).
Does it appear \(\sigma_{Y|X}\) depends upon X? Not much evidence for this.
Do you see obvious outliers? Observation 3 is an outlier in the \(x\) direction, and therefore possibly a high leverage point and influential on the model fit.
Is the normality assumption reasonable? There appears to be some skewness, but with so few points normality may be reasonable.
Is there a striking pattern in residuals vs. order of the data? No striking pattern.
Check the Cook’s \(D\) values. We anticipated that the \(3^{\textrm{rd}}\) observation is affecting the fit by a lot more than any other values. The \(D\)-value is much larger than 1. Note that the residual is not large for this value.
Omit problem observations from the analysis and see if any conclusions change substantially. Let us refit the equation without observation 3 to see if anything changes drastically. I will use the weighted least squares approach discussed earlier on this example. Define a variable wt that is 1 for all observations except obs. 3, and make it 0 for that one.
# wt = 1 for all except obs 3 where wt = 0dat_thyroid <- dat_thyroid |>mutate(wt =ifelse(id ==3, 0, 1) )dat_thyroid$wt
[1] 1 1 0 1 1 1 1 1
What changes by deleting case 3? The fitted line gets steeper
(slope changes from $-1.30$ to $-2.19$), adjusted $R^2$ gets larger
(up to 58% from 53%), and $S$ changes from 11.7 to 10.6.
Because the Weight values are much less spread out,
$SE(\hat{\beta_1})$ becomes quite a bit larger (to 0.714, up from 0.436)
and we lose a degree of freedom for MS Error (which will
penalize us on tests and CIs). Just about any quantitative
statement we would want to make using CIs would be about the same
either way since CIs will overlap a great deal, and our
qualitative interpretations are unchanged (Blood Loss drops with
Weight). Unless something shows up in the plots, I don't see any
very important changes here.
lm_blood_wt_no3 <-lm(blood_loss ~ weight, data = dat_thyroid, weights = wt)# use summary() to get t-tests of parameters (slope, intercept)summary(lm_blood_wt_no3)
# exclude obs 3dat_thyroid_no3 <- dat_thyroid |>filter( wt ==1 )# ggplot: Plot the data with linear regression fit and confidence bandslibrary(ggplot2)p <-ggplot(dat_thyroid, aes(x = weight, y = blood_loss, label = id))p <- p +geom_point()# plot labels next to pointsp <- p +geom_text(hjust =0.5, vjust =-0.5)# plot regression line and confidence bandp <- p +geom_smooth(method = lm, color ="blue", se =FALSE, linetype ="solid")p <- p +geom_smooth(data = dat_thyroid_no3, method = lm, color ="red", se =FALSE, linetype ="dashed", fullrange =TRUE)p <- p +labs(caption ="Blue solid is full data, Red dashed is without obs 3.")print(p)
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Nothing very striking shows up in the residual plots, and no
Cook's $D$ values are very large among the remaining observations.
Error in qr.resid(xqr, yt) :
'qr' and 'y' must have the same number of rows
How much difference is there in a practical sense? Examine the 95% prediction interval for a new observation at Weight = 50kg. Previously we saw that interval based on all 8 observations was from 457.1 to 517.8 ml of Blood Loss. Based on just the 7 observations the prediction interval is 451.6 to 512.4 ml. There really is no practical difference here.
# CI for the mean and PI for a new observation at weight=50predict(lm_blood_wt , data.frame(weight=50), interval ="prediction")
Warning in predict.lm(lm_blood_wt_no3, data.frame(weight = 50), interval = "prediction"): Assuming constant prediction variance even though model fit is weighted
fit lwr upr
1 481.9939 451.5782 512.4096
Therefore, while obs. 3 was potentially influential, whether the value is included or not makes very little difference in the model fit or relationship between Weight and BloodLoss.
Gesell data
These data are from a UCLA study of cyanotic heart disease in children. The predictor is the age of the child in months at first word and the response variable is the Gesell adaptive score, for each of 21 children.
'data.frame': 21 obs. of 3 variables:
$ id : int 1 2 3 4 5 6 7 8 9 10 ...
$ age : int 15 26 10 9 15 20 18 11 8 20 ...
$ score: int 95 71 83 91 102 87 93 100 104 94 ...
Plot Score versus Age. Comment on the relationship between Score and Age.
# ggplot: Plot the data with linear regression fit and confidence bandslibrary(ggplot2)p <-ggplot(dat_gesell, aes(x = age, y = score, label = id))p <- p +theme_bw()p <- p +geom_point()# plot labels next to pointsp <- p +geom_text(hjust =0.5, vjust =-0.5)# plot regression line and confidence band#p <- p + geom_smooth(method = lm)p <- p +geom_smooth(method = lm, color ="blue", se =FALSE, linetype ="solid")p <- p +geom_smooth(data = dat_gesell |>filter(!(id ==18)), method = lm, color ="red", se =FALSE, linetype ="dashed", fullrange =TRUE)p <- p +geom_smooth(data = dat_gesell |>filter(!(id ==19)), method = lm, color ="purple", se =FALSE, linetype ="dotted", fullrange =TRUE)p <- p +geom_smooth(data = dat_gesell |>filter(!(id %in%c(18, 19))), method = lm, color ="orange", se =FALSE, linetype ="solid", fullrange =TRUE)p <- p +labs(caption ="Blue solid is full data, Red dashed is without obs 18, Green dotted is without obs 19.")print(p)
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
There are no obvious transformations to try here.
Fit a simple linear regression model. Provide an equation for the LS line. Does age at first word appear to be an “important predictor” of Gesell adaptive score? (i.e., is the estimated slope significantly different from zero?)
lm_score_age <-lm(score ~ age, data = dat_gesell)# use summary() to get t-tests of parameters (slope, intercept)summary(lm_score_age)
Call:
lm(formula = score ~ age, data = dat_gesell)
Residuals:
Min 1Q Median 3Q Max
-15.604 -8.731 1.396 4.523 30.285
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 109.8738 5.0678 21.681 7.31e-15 ***
age -1.1270 0.3102 -3.633 0.00177 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11.02 on 19 degrees of freedom
Multiple R-squared: 0.41, Adjusted R-squared: 0.3789
F-statistic: 13.2 on 1 and 19 DF, p-value: 0.001769
Do these plots suggest any inadequacies with the model?
# ggplot: Plot the data with linear regression fit and confidence bandslibrary(ggplot2)p <-ggplot(dat_gesell |>filter(age <=20), aes(x = age, y = score, label = id))p <- p +theme_bw()p <- p +geom_point()# plot labels next to pointsp <- p +geom_text(hjust =0.5, vjust =-0.5)# plot regression line and confidence band#p <- p + geom_smooth(method = lm)p <- p +geom_smooth(method = lm, color ="blue", se =FALSE, linetype ="solid")print(p)
`geom_smooth()` using formula = 'y ~ x'
Warning: The following aesthetics were dropped during statistical transformation: label
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
