hospital_data <- read.csv("C:/Users/ajste/Downloads/week 6 data.csv")

1.

Using R, conduct correlation analysis (between the two variables) and interpret.

# Correlation test
cor.test(hospital_data$RVUs, hospital_data$Expenditures)
## 
##  Pearson's product-moment correlation
## 
## data:  hospital_data$RVUs and hospital_data$Expenditures
## t = 46.449, df = 382, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9051402 0.9355065
## sample estimates:
##       cor 
## 0.9217239
# Scatter-plot graph of correlation
x <- hospital_data$RVUs / 10^6
y <- hospital_data$Expenditures / 10^9

plot(x, y,
     main = "Correlation between RVUs and Expenditures",
     xlab = "RVUs (x 10^6)",
     ylab = "Expenditures (x $10^9 USD)",
     col = "blue",
     pch = 19,
     cex = 1.5)

abline(lm(y ~ x), col = "red", lwd = 2)

r_label  <- bquote(italic(r) == .(0.9217239))

legend("bottomright", 
       legend = c(as.expression(r_label)), 
       bty = "n",      
       text.col = "black", 
       cex = 1.2)     

For this correlation analysis, the null hypothesis states that there is no correlation between RVUs and expenditures while the alternative hypothesis states that there is a significant correlation between RVUs and expenditures. Running a t-test results in a t-value of 46.449 which has a very small p-value of < 2.2\(e^{-16}\), which indicates that the relationship is statistically significant. Given a correlation coefficient of 0.9217239, we can say that there is a very strong positive correlation between a hospital’s RVUs and their expenditures. A strong correlation between RVUs and expenditures makes sense as greater numbers or intensities of service would also indicate a greater cost required to perform those services, however this does not necessarily mean that increasing RVUs is the direct causation of increasing expenditures.

2.

Then fit a linear model with Expenditure as the dependent variable (Y) and RVUs as the independant (X) variable. Interpret the results (Interpreting regression coefficients in particular) and check whether the Gauss Markov Assumptions / linear regression assumptions hold or not (by conducting residual plot analysis and explaining your results in your own words).

# Fit linear regression model
linear_model <- lm(Expenditures ~ RVUs, data = hospital_data)

# Regression results
summary(linear_model)
## 
## Call:
## lm(formula = Expenditures ~ RVUs, data = hospital_data)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -185723026  -14097620    2813431   11919781  642218316 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.785e+06  4.413e+06  -0.858    0.392    
## RVUs         2.351e+02  5.061e+00  46.449   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 67350000 on 382 degrees of freedom
## Multiple R-squared:  0.8496, Adjusted R-squared:  0.8492 
## F-statistic:  2157 on 1 and 382 DF,  p-value: < 2.2e-16
# Coefficients for regression equation
coef(linear_model)
##   (Intercept)          RVUs 
## -3785072.1584      235.0719
# Residual diagnostic plots
par(mfrow = c(2, 2))
plot(linear_model)

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
# Breusch-Pagan test (constant variance)
bptest(linear_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  linear_model
## BP = 62.072, df = 1, p-value = 3.311e-15
# Shapiro-Wilks test (Normality)
shapiro.test(residuals(linear_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(linear_model)
## W = 0.62591, p-value < 2.2e-16
# Durbin-Watson test (autocorrelation)
dwtest(linear_model)
## 
##  Durbin-Watson test
## 
## data:  linear_model
## DW = 1.9515, p-value = 0.3097
## alternative hypothesis: true autocorrelation is greater than 0

The estimated regression equation for the relationship between RVUs and expenditures is Expenditures = -3.7851\(*10^{6}\) + 235.0719(RVUs). A coefficient of 235.0719 for RVUs means that a one-unit increase in RVUs is associated with an increase of approximately $235.07 in expenditures. A p-value of < 2\(e^{-16}\) indicates that the RVU coefficient is statistically significant. The intercept, -3.7851\(*10^{6}\), has a p-value of 0.392, which indicates that this coefficient is not statistically significant, which makes sense because a hospital with zero RVUs is unrealistic. A \(R^{2}\) value of 0.8492 indicates that approximately 85% of the variation in expenditures can be explained by the variation in RVUs.

The residuals-versus-fitted plot shows that the majority of points are clustered around a flat, horizontal line at zero, which would indicate a perfectly linear relationship. However, the spread of the residuals becomes much larger as the fitted expenditure values increase. This suggests that an assumption of a perfectly linear relationship is questionable and that a simple linear model does not completely explain the structure seen in the data. The Q-Q residuals plot shows the upper tail distancing from the reference line with several large positive residuals. This indicates the distribution of residuals is strongly right-skewed, violating the assumption of normality. The scale-location plot shows an increase in the spread of residuals as the fitted expenditure values increase, indicating heteroscedasticity in the linear model. The residuals-versus-leverage plot shows that while most points are clustered around zero, there are several observations with high Cook’s distance values which have a large influence on the model’s fitted values. Overall, the linear regression assumptions are not held. The non-constant variance and non-normality of the residuals is also backed up by the results of the Breusch-Pagan and Shapiro-Wilk tests.

3.

Then fit a linear model of ln(Expenditures)~RVUs. Mathematically speaking, the logarithm function tends to squeeze together the larger values in your data set and stretches out the smaller values. (If you are wondering why do a log transformation, see the first two charts here that shows how log reduces skewness to help meet normality of X assumption, with the caveat being that X is somewhat normally distributed to begin with in order for the transformation to reduce / remove skewness). This transformation is routine in Economics or Finance forecasting to stabilize the variance of a time series (GDP, stock prices,…).

# Fit log-transformed model
log_model <- lm(log(Expenditures) ~ RVUs, data = hospital_data)

# Log-transformed regression results
summary(log_model)
## 
## Call:
## lm(formula = log(Expenditures) ~ RVUs, data = hospital_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.59439 -0.29504  0.06135  0.35333  1.20871 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.730e+01  3.325e-02  520.11   <2e-16 ***
## RVUs        1.349e-06  3.814e-08   35.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5076 on 382 degrees of freedom
## Multiple R-squared:  0.7661, Adjusted R-squared:  0.7655 
## F-statistic:  1251 on 1 and 382 DF,  p-value: < 2.2e-16
# Coefficients for log-transformed regression equation
coef(log_model)
##  (Intercept)         RVUs 
## 1.729584e+01 1.349109e-06
# Residual diagnostic plots
par(mfrow = c(2, 2))
plot(log_model)

# Breusch-Pagan test (constant variance)
bptest(log_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  log_model
## BP = 25.485, df = 1, p-value = 4.458e-07
# Shapiro-Wilks test (Normality)
shapiro.test(residuals(log_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(log_model)
## W = 0.98687, p-value = 0.001513
# Durbin-Watson test (autocorrelation)
dwtest(log_model)
## 
##  Durbin-Watson test
## 
## data:  log_model
## DW = 2.1114, p-value = 0.8587
## alternative hypothesis: true autocorrelation is greater than 0

How did this log transformation affect the Gauss Markov Assumptions (sure, the residual analysis diagnostic charts will change but what is your takeaway - which assumptions are better met or are some assumptions not met now)?

The log-transformed model exhibits much greater normality than the standard linear regression model. Unlike in the linear model, the Q-Q residuals plot does not have any tails that deviate from the reference line. While less skewed than the linear model, the log-transformed is still not normal enough to satisfy the assumption of normality, with the Shapiro-Wilk test remaining statistically significant. The log-transformed model also exhibits greater homoscedasticity than the linear model. Unlike in the linear model, the residuals-versus-leverage plot does not show the same fanning outward from the reference line as the fitted expenditure values increase. However, the log-transformed model still doesn’t meet the assumption of constant variance, with a statistically significant Breusch-Pagan test. The maximum Cook’s distance also decreases substantially between the linear and log-transformed models, meaning that the influence of large observations is reduced.

Are you happy with this functional form capturing the relationship between Y and X or would like to keep some different functional form (EG - ln(Expenditures)~ln(RVUs) or ln(Expenditures)~RVUs + RVUs^2)? Why (4 lines maximum)? In other words, would you put all your money to create hospital expansion plans (hiring more doctors and nurses, opening more rooms,…) based on the specific relationship you find between expenditure and RVUs (or you think it could blow up in your face)? You can try some other transformations for dealing with positively skewed data like root to the power n or even reciprocal.

# Fit double log-transformed model
log_log_model <- lm(log(Expenditures) ~ log(RVUs), data = hospital_data)

# Double log-transformed regression results
summary(log_log_model)
## 
## Call:
## lm(formula = log(Expenditures) ~ log(RVUs), data = hospital_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74657 -0.19864 -0.02431  0.18642  0.93551 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.91487    0.16621   41.60   <2e-16 ***
## log(RVUs)    0.88444    0.01317   67.17   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2932 on 382 degrees of freedom
## Multiple R-squared:  0.9219, Adjusted R-squared:  0.9217 
## F-statistic:  4512 on 1 and 382 DF,  p-value: < 2.2e-16
# Coefficients for double log-transformed regression equation
coef(log_log_model)
## (Intercept)   log(RVUs) 
##   6.9148650   0.8844391
# Residual diagnostic plots
par(mfrow = c(2, 2))
plot(log_log_model)

# Breusch-Pagan test (constant variance)
bptest(log_log_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  log_log_model
## BP = 0.74786, df = 1, p-value = 0.3872
# Shapiro-Wilk test (normality)
shapiro.test(residuals(log_log_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(log_log_model)
## W = 0.98887, p-value = 0.005059
# Durbin-Watson test (autocorrelation)
dwtest(log_log_model)
## 
##  Durbin-Watson test
## 
## data:  log_log_model
## DW = 1.7349, p-value = 0.00436
## alternative hypothesis: true autocorrelation is greater than 0

A ln(Expenditures) ~ ln(RVUs) model, I believe, is a much better functional form of this relationship. This model eliminates evidence of heteroscedasticity, which is supported by a Breusch-Pagan test that is not statistically significant (P = 0.3872). This model also has a much higher \(R^{2}\) value (0.9217) than the ln(Expenditures) ~ RVUs model (0.7655). However, the residuals of this model are still not perfectly normal, so I would still be hesitant about putting all my money on hospital expansion plans.