Newey West Standard error adjustment in GLM

Question

I am trying to solve this problem since weeks now and I am a little desperate. I am new to R so please be patient with me. Thus, help is highly appreciated!

Here is the thing:

I am running a regression discontinuity in time analysis.
I have a gamma distributed outcome variable.
I have autocorrelated residuals with a lag of 2
I want to plot my regression line with Newey West SE adjusted CI and PI

--> See my output: plot

I have several questions:

Does this code make sense? Is the Newey West adjustment correct? Are there simpler approaches?
Do you know how I can smooth the CI?
The CI appears to wide, do you have an idea, why?

Thanks in advance, and sorry, if something is wrongly coded.plot

This is my code:

#Packages
------------------------------------------------------------------
library(jtools)
library(sandwich)

#Calculation
------------------------------------------------------------------
# GLM
model <- glm(outcome ~ time * x1, data = dat, family = Gamma(link = "log"))

# Compute Newey-West standard errors with the specified lag order = 2
vcovNW_model <- NeweyWest(model, lag = 2)

# Perform coefficient test with Newey-West standard errors
NW_model <- coeftest(model, vcov. = vcovNW_model)

# Compute confidence intervals
CI_NW_model <- exp(confint(NW_model))

# Calculate predicted values and standard errors
PRED_model <- predict(model, newdata = dat, type = "response", se.fit = TRUE)

# Extract predictions and standard errors
PREDICTIONS_model <- PRED_model$fit
SE_model <- exp(NW_model[, "Std. Error"])

# Calculate lower and upper bounds of the confidence interval
LOW_CI <- PREDICTIONS_model - 1.96 * SE_model
UP_CI  <- PREDICTIONS_model + 1.96 * SE_model

# Calculate prediction interval
PI_model <- 1.96 * sqrt(PRED_model$se.fit^2 + SE_model^2)

#Plot
------------------------------------------------------------------
# Plot axes
plot(0, type = "n",
     xlim = c(-19, 29),
     ylim = c(0, 15),
     xlab = "x1",
     ylab = "Outcome",
     xaxt = "n")
axis(side = 1, at = seq(-19, 29, by = 5))

# Plot prediction interval as a gray shaded area
polygon(
  c(dat$x1, rev(dat$x1)),
  c(PREDICTIONS_model - PI_model, rev(PREDICTIONS_model + PI_model)),
col =  rgb(0, 153, 102, maxColorValue = 255, alpha = 50), border = NA, alpha = 0.3
)

# Plot confidence interval as a shaded area
polygon(
  c(lolli_SR$x1, rev(lolli_SR$x1)),
  c(LOW_CI, rev(UP_CI)),
  col =  rgb(0, 153, 102, maxColorValue = 255, alpha = 100), border = NA, alpha = 0.3
)

# Plot regression line
lines(lolli_SR$x1, PREDICTIONS_model, col = "#009966", lwd = 2)


# Plot data points 
points(lolli_SR$x1, lolli_SR$outcome, col = "black", pch = 21, bg = "grey", cex = 0.75)

Newey West Standard error adjustment in GLM

0 Answers0