147

library(AER)

## Loading required package: car

## Loading required package: carData

## Loading required package: lmtest

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: sandwich

## Loading required package: survival

data("DoctorVisits")
d<-DoctorVisits
multiple_regression_analysis <- function(data, dependent_var, independent_vars) {
  formula <- paste(dependent_var, "~", paste(independent_vars, collapse = "+"))
  model <- lm(formula, data = data)
  prediction <- predict(model, newdata = data)

  bp_test <- bptest(model)
  if (is.na(bp_test$p.value)) {
    cat("Breusch-Pagan test failed. Check your data for heteroskedasticity.")
  } else {
    p_value <- bp_test$p.value
    predictions <- predict(model, newdata = data)
  
    return(list(model = model, p_value = p_value, predictions = predictions))
  }
}

Y <- "visits"
x <- c("age", "income", "illness", "reduced", "health", "private", "freepoor", "freerepat", "nchronic", "lchronic")

result <- multiple_regression_analysis(d, Y, x)

summary(result$model)

## 
## Call:
## lm(formula = formula, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1321 -0.2555 -0.1440 -0.0487  7.0513 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.050523   0.034605   1.460 0.144363    
## age           0.158125   0.066524   2.377 0.017492 *  
## income       -0.066550   0.030394  -2.190 0.028597 *  
## illness       0.060295   0.008336   7.233  5.4e-13 ***
## reduced       0.103182   0.003657  28.218  < 2e-16 ***
## health        0.017212   0.005179   3.323 0.000896 ***
## privateyes    0.041578   0.024547   1.694 0.090367 .  
## freepooryes  -0.105316   0.052446  -2.008 0.044687 *  
## freerepatyes  0.038542   0.037905   1.017 0.309286    
## nchronicyes   0.006572   0.023683   0.277 0.781426    
## lchronicyes   0.041868   0.035532   1.178 0.238724    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7139 on 5179 degrees of freedom
## Multiple R-squared:  0.2014, Adjusted R-squared:  0.1999 
## F-statistic: 130.6 on 10 and 5179 DF,  p-value: < 2.2e-16

result$p_value

##            BP 
## 1.757476e-146

head(result$predictions)

##         1         2         3         4         5         6 
## 0.5757780 0.3760682 0.2015555 0.1308785 0.7109039 0.6234081