Estimate confidence intervals from a model with clustered standard errors

Question

I am trying to find a way to estimate the predicted values of y with confidence intervals for specific values of x in an OLS regression. My model includes an interaction term and I use clustered standard errors and weights in my model.

A similar question was asked and answer previously, I thought it could be a good starting point:

robust standard errors in ggplot2

The problem is that, the solution offered here does not work when there are interaction terms OR weights in the model. It does produces an outcome when there are both weights and interaction terms. I found this confusing but I am relatively new to R and I could not understand the source of the problem.

In the second and third examples (lm2 & lm3) I get "Error in X %*% V : non-conformable arguments". My best guess for the source of the error in the third case is that model.frame(lm3) does not include interaction terms. But I don’t know whether I am in the right track and could not find a way to fix it. Besides it's not clear to me how I can set x1 to a specific value in this example. Can someone help me revise to code above or offer an alternative way to get my predicted standard errors when x is set to a specific value ?

df <- data.frame(x1 = rnorm(100), x2 = rnorm(100), w1 = runif(100,0.1,2),y = rnorm(100), group = as.factor(sample(1:10, 100, replace=T))) 
        lm1 <- lm(y ~ x1+x2, data = df)
        lm2 <- lm(y ~ x1+x2, data = df, weight=w1)
        lm3 <- lm(y ~ x1*x2, data = df)
        lm4 <- lm(y ~ x1*x2, data = df, weight=w1)

        getvcov <- function(fm,dfcw,cluster) {
          library(sandwich);library(lmtest)
          M <- length(unique(cluster))   
          N <- length(cluster)           
          K <- fm$rank                        
          dfc <- (M/(M-1))*((N-1)/(N-K))  
          uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
          dfc*sandwich(fm, meat=crossprod(uj)/N)*dfcw
        }

        V <- getvcov(lm1,1,df$group)
        X <- as.matrix(model.frame(lm1))
        se <- predict(lm1,se=TRUE)$se.fit
        se_robust1 <- sqrt(diag(X %*% V %*% t(X)))

        V <- getvcov(lm2,1,df$group)
        X <- as.matrix(model.frame(lm2))
        se <- predict(lm2,se=TRUE)$se.fit
        se_robust2 <- sqrt(diag(X %*% V %*% t(X)))

        V <- getvcov(lm3,1,df$group)
        X <- as.matrix(model.frame(lm3))
        se <- predict(lm3,se=TRUE)$se.fit
        se_robust2 <- sqrt(diag(X %*% V %*% t(X)))

        V <- getvcov(lm4,1,df$group)
        X <- as.matrix(model.frame(lm4))
        se <- predict(lm4,se=TRUE)$se.fit
        se_robust4 <- sqrt(diag(X %*% V %*% t(X)))

Estimate confidence intervals from a model with clustered standard errors

0 Answers0