Introduction

In this project, I’m going to analyze the relationships between patient satisfaction and several patient-related factors using linear regression.

Loading Data

The data for this project comes from Applied Linear Statistical Models, Fifth Edition, by Michael H. Kutner, Christopher J. Nachtsheim, John Neter, and William Li.

satisfactiondata <- read.table("http://www.cnachtsheim-text.csom.umn.edu/Kutner/Chapter%20%206%20Data%20Sets/CH06PR15.txt", sep ="" , header = FALSE)

colnames(satisfactiondata) <- c("Satisfaction", "Age", "IllnessSeverity", "Anxiety")

Dataset Context:

A hospital administrator wished to study the relation between patient satisfaction (Y) and patient’s age (X I, in years), severity of illness (X2, an index), and anxiety level (X3, an index). The administrator randomly selected 46 patients and collected the data presented below, where larger values of Y, X2 , and X3 are, respectively, associated with more satisfaction, increased severity of illness, and more anxiety.

Research Questions

I’m going to use research questions to guide my analysis.

Question 1. How do age, severity of illness, and anxiety level collectively predict patient satisfaction?

Question 2. Which factors are significantly associated with patient satisfaction after controlling for the other variables?

Question 3. Are there interaction effects between predictors that influence patient satisfaction?

Data Analysis

Loading Packages

library(dplyr)
library(car)
library(lmtest)

Data Exploration

First, I’d like to get a glimpse of what the dataset looks like.

glimpse(satisfactiondata)
## Rows: 46
## Columns: 4
## $ Satisfaction    <int> 48, 57, 66, 70, 89, 36, 46, 54, 26, 77, 89, 67, 47, 51…
## $ Age             <int> 50, 36, 40, 41, 28, 49, 42, 45, 52, 29, 29, 43, 38, 34…
## $ IllnessSeverity <int> 51, 46, 48, 44, 43, 54, 50, 48, 62, 50, 48, 53, 55, 51…
## $ Anxiety         <dbl> 2.3, 2.3, 2.2, 1.8, 1.8, 2.9, 2.2, 2.4, 2.9, 2.1, 2.4,…

To further explore the data, I’m going to create stem-and-leaf plots for each predictor variable to examine their distributions.

cat("Stem-and-Leaf Plot for Age\n")
## Stem-and-Leaf Plot for Age
stem(satisfactiondata$Age)
## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   2 | 23
##   2 | 58899999
##   3 | 012233344
##   3 | 66678
##   4 | 0012233344
##   4 | 557779
##   5 | 0233
##   5 | 55
cat("\nStem-and-Leaf Plot for Illness Severity\n")
## 
## Stem-and-Leaf Plot for Illness Severity
stem(satisfactiondata$IllnessSeverity)
## 
##   The decimal point is at the |
## 
##   40 | 0
##   42 | 00
##   44 | 0
##   46 | 00000
##   48 | 0000000000
##   50 | 000000000000
##   52 | 000000
##   54 | 0000
##   56 | 00
##   58 | 0
##   60 | 0
##   62 | 0
cat("\nStem-and-Leaf Plot for Anxiety\n")
## 
## Stem-and-Leaf Plot for Anxiety
stem(satisfactiondata$Anxiety)
## 
##   The decimal point is 1 digit(s) to the left of the |
## 
##   18 | 000000
##   20 | 00000000
##   22 | 0000000000000
##   24 | 00000000000
##   26 | 0000
##   28 | 0000

The plot for age appears roughly bimodal, with peaks around the early 30s and early 40s. The plots for illness severity and anxiety appear approximately symmetric.

Model Fitting

I’m going to create a multiple linear regression model with patient satisfaction as the dependent variable and age, anxiety, and severity of illness as predictor variables. This model will allow me to examine the relationships between patient satisfaction and each predictor while controlling for the others.

model <- lm(Satisfaction ~ Age + IllnessSeverity + Anxiety, data = satisfactiondata) 

summary(model)
## 
## Call:
## lm(formula = Satisfaction ~ Age + IllnessSeverity + Anxiety, 
##     data = satisfactiondata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.3524  -6.4230   0.5196   8.3715  17.1601 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     158.4913    18.1259   8.744 5.26e-11 ***
## Age              -1.1416     0.2148  -5.315 3.81e-06 ***
## IllnessSeverity  -0.4420     0.4920  -0.898   0.3741    
## Anxiety         -13.4702     7.0997  -1.897   0.0647 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.06 on 42 degrees of freedom
## Multiple R-squared:  0.6822, Adjusted R-squared:  0.6595 
## F-statistic: 30.05 on 3 and 42 DF,  p-value: 1.542e-10

This output shows that at alpha 0.05, the coefficient for Age is significant, with a p-value of p < 0.0001. The coefficient for Anxiety is approaching significance, with a p-value of 0.0647. The coefficient for IllnessSeverity is not significant, with a p-value of 0.3741.

The adjusted R-squared of 0.6595 indicates a good model fit, with the predictors explaining approximately 66% of the variability in patient satisfaction scores.

The F-statistic (30.05, p < 0.0001) also indicates that the overall model is highly significant and provides a strong fit to the data.

The estimated coefficient for Age shows that, holding illness severity and anxiety constant, each additional year of age is associated with a decrease of approximately 1.14 units in patient satisfaction. Similarly, the coefficient for Anxiety suggests that, for each one-point increase in anxiety, patient satisfaction decreases by roughly 13.5 units, though this effect is only approaching significance.

Testing Linear Model Assumptions

Now, I need to check the assumptions of the linear model to ensure this model is a good fit.

The assumptions of a multiple linear regression model are as follows:

1. Linearity

2. Homoskedasticity

3. Independence of Error

4. Normality

5. Independence of the Independent Variables (multicollinearity)

Testing Linearity and Homoskedasticity

I’m going to use a residual plot to test the linearity and homoskedasticity assumptions.

plot(model, which = 1)

The residual plot shows no violations of the model assumptions. The red line is nearly horizontal, indicating the linearity assumption is met The residuals indicate roughly constant variance, indicating that the homoskedasticity assumption is also met.

I’m also going to use a Breusch-Pagan test to test for homoskedasticity to confirm.

bptest(model)
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 2.5583, df = 3, p-value = 0.4648

The Breusch-Pagan test identifies a p-value of 0.4648 which is not significant. I’ll fail to reject the null hypothesis, indicating that there is not sufficient evidence of heteroskedasticity in the model. The homoskedasticity assumption appears to be met.

Independence of Errors

The independence of errors assumption is assumed to be met because the patients were independently sampled, and the data does not involve time-series or repeated measurements. If this data was time-ordered, I would use a Durbin-Watson test to evaluate the independence of errors.

Testing Normality

Now I’m going to use a normal probability plot and a histogram to test the normality assumption.

plot(model, which=2)

The normal probability plot looks nearly normal with slight deviation in the upper and lower tails.

hist(residuals(model),
     breaks = 10)

This histogram appears nearly normal with a slight left skew.

These graphs indicate the normality assumption is met.

Testing Independence of the Independent Variables

Finally, I’m going to test the Independence of the Independent Variables assumption using the variance inflated factor.

vif(model)
##             Age IllnessSeverity         Anxiety 
##        1.632296        2.003235        2.009062

All VIF values appear to be below 5, indicating there is no concerning evidence of multicollinearity among the predictor variables.

Assumptions Conclusion

The Linearity, Homoskedasticity, Independence of Error, Normality, and Independence of the Independent Variables assumptions all appear to be met, indicating that the linear model is appropriate for this dataset.

Interaction Terms

model <- lm(Satisfaction ~ Age + IllnessSeverity + Anxiety + Age * IllnessSeverity, data = satisfactiondata) 

summary(model)
## 
## Call:
## lm(formula = Satisfaction ~ Age + IllnessSeverity + Anxiety + 
##     Age * IllnessSeverity, data = satisfactiondata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.2971  -5.9951   0.3513   8.4263  17.2372 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)  
## (Intercept)         135.70057   79.72133   1.702   0.0963 .
## Age                  -0.55257    2.01699  -0.274   0.7855  
## IllnessSeverity       0.01017    1.61770   0.006   0.9950  
## Anxiety             -13.44757    7.17858  -1.873   0.0682 .
## Age:IllnessSeverity  -0.01159    0.03944  -0.294   0.7704  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.17 on 41 degrees of freedom
## Multiple R-squared:  0.6829, Adjusted R-squared:  0.6519 
## F-statistic: 22.07 on 4 and 41 DF,  p-value: 8.946e-10

The interaction term between Age and IllnessSeverity is not significant (p = 0.7704).

model <- lm(Satisfaction ~ Age + IllnessSeverity + Anxiety + Age * Anxiety, data = satisfactiondata) 

summary(model)
## 
## Call:
## lm(formula = Satisfaction ~ Age + IllnessSeverity + Anxiety + 
##     Age * Anxiety, data = satisfactiondata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.8727  -6.7365   0.8366   9.0390  15.7261 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     104.7373    55.8137   1.877   0.0677 .
## Age               0.2245     1.3587   0.165   0.8696  
## IllnessSeverity  -0.4013     0.4934  -0.813   0.4206  
## Anxiety           9.5311    23.6778   0.403   0.6894  
## Age:Anxiety      -0.5975     0.5868  -1.018   0.3145  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.05 on 41 degrees of freedom
## Multiple R-squared:   0.69,  Adjusted R-squared:  0.6598 
## F-statistic: 22.82 on 4 and 41 DF,  p-value: 5.653e-10

The interaction term between Age and Anxiety is not significant (p = 0.3145).

model <- lm(Satisfaction ~ Age + IllnessSeverity + Anxiety + IllnessSeverity * Anxiety, data = satisfactiondata) 

summary(model)
## 
## Call:
## lm(formula = Satisfaction ~ Age + IllnessSeverity + Anxiety + 
##     IllnessSeverity * Anxiety, data = satisfactiondata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.4007  -5.8432   0.4612   8.5402  17.0120 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             140.14938  110.40721   1.269    0.211    
## Age                      -1.14052    0.21742  -5.246 5.07e-06 ***
## IllnessSeverity          -0.07586    2.22961  -0.034    0.973    
## Anxiety                  -5.55837   47.50887  -0.117    0.907    
## IllnessSeverity:Anxiety  -0.15716    0.93287  -0.168    0.867    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.18 on 41 degrees of freedom
## Multiple R-squared:  0.6824, Adjusted R-squared:  0.6514 
## F-statistic: 22.02 on 4 and 41 DF,  p-value: 9.203e-10

The interaction term between Age and Anxiety is not significant (p = 0.867).

Answering Research Questions

Question 1. How do age, severity of illness, and anxiety level collectively predict patient satisfaction?

I used a multiple linear regression model to examine the effects of the predictor variables on patient satisfaction. The F-statistic indicated that the model was highly significant (30.05, p < 0.0001). The adjusted r-squared of 0.6595 indicates that approximately 66% of the variability in patient satisfaction can be explained by the three predictors.

Question 2. Which factors are significantly associated with patient satisfaction after controlling for the other variables?

When controlling for other predictors in the model,

  • Age significantly predicted patient satisfaction (p < 0.0001). Older patients tended to report lower satisfaction, with each additional year of age being associated with an approximate 1.14 unit decrease in satisfaction.
  • Anxiety showed a negative relationship with satisfaction, with higher anxiety associated with lower satisfaction, but this effect was only approaching significance (p = 0.0647). The approximate decrease in satisfaction is about 13.5 units per one-point increase in anxiety.
  • Illness severity was not significantly associated with patient satisfaction after controlling for age and anxiety (p = 0.3741).

Overall, age was the strongest and only significant predictor, anxiety may have had a moderate effect, and illness severity did not meaningfully contribute to predicting patient satisfaction.

Question 3. Are there interaction effects between predictors that influence patient satisfaction?

To examine the interactions between predictor variables, I tested three interaction terms: Age x IllnessSeverity, Age x Anxiety, and IllnessSeverity x Age. None of these interaction variables significantly influenced patient satisfaction, so there were no interaction effects that influenced patient satisfaction.

Conclusion

In conclusion, patient satisfaction was significantly influenced by age, with older patients reporting lower satisfaction. Anxiety may have a slight negative effect, while illness severity does not appear to impact satisfaction when controlling for the other factors. No significant interaction effects were found among the predictors, suggesting that the relationships between age, anxiety, and illness severity with satisfaction are largely independent. Overall, the multiple linear regression model explained a substantial portion of the variability in patient satisfaction and provides insight into which patient characteristics are most relevant for predicting satisfaction levels.