3

我正在尝试找到一种方法来估计 y 的预测值，并使用 OLS 回归中 x 的特定值的置信区间。我的模型包含一个交互项，并且我在模型中使用了聚类标准误差和权重。

之前提出并回答了一个类似的问题，我认为这可能是一个很好的起点：

ggplot2中的稳健标准误差

问题是，当模型中有交互项或权重时，此处提供的解决方案不起作用。当同时存在权重和交互项时，它确实会产生结果。我发现这令人困惑，但我对 R 比较陌生，我无法理解问题的根源。

在第二个和第三个示例（lm2 和 lm3）中，我得到“X %*% V 中的错误：不符合要求的参数”。在第三种情况下，我对错误来源的最佳猜测是 model.frame(lm3) 不包括交互项。但我不知道我是否走在正确的轨道上，也找不到解决方法。此外，我不清楚如何在此示例中将 x1 设置为特定值。当 x 设置为特定值时，有人可以帮我修改上面的代码或提供另一种方法来获得我的预测标准错误吗？

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)))
4

0 回答 0