Untitled

#lấy dữ liệu và chạy hồi quy

library(AER)
data("DoctorVisits")
d<- DoctorVisits
X<- d$age
Y<- d$income
regression_heteroskedasticity_forecast <- function(data, y_var, x_var) {

  model <- lm(formula = paste(y_var, "~", x_var), data = data)

  bp_test <- bptest(model)
  p_value <- bp_test$p.value

  prediction <- predict(model, newdata = data)

  return(list(model = model, p_value = p_value, prediction = prediction))
}

result <- regression_heteroskedasticity_forecast(data.frame(X, Y), "Y", "X")

summary(result$model)

## 
## Call:
## lm(formula = paste(y_var, "~", x_var), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.68883 -0.22885 -0.09976  0.22582  1.06999 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.78161    0.01096   71.35   <2e-16 ***
## X           -0.48833    0.02407  -20.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3551 on 5188 degrees of freedom
## Multiple R-squared:  0.07348,    Adjusted R-squared:  0.0733 
## F-statistic: 411.5 on 1 and 5188 DF,  p-value: < 2.2e-16

result$p_value

##           BP 
## 8.864352e-13

head(result$prediction)

##        1        2        3        4        5        6 
## 0.688827 0.688827 0.688827 0.688827 0.688827 0.688827

```