r - R 中看似不相关的回归与估算数据池结果

Question

我正在尝试使用 R 中的 systemfit 包完成看似不相关的回归 (SUR)。但是，使用多重插补数据(使用 mouse 包)完成这些分析并不简单。

在谷歌上搜索这个问题时，我发现有一篇关于相同问题的已删除帖子，该帖子似乎使用了以下示例(归功于发布者，进行了少量编辑)

library(systemfit)
library(mice)
nhanes2

r1 <- bmi ~ hyp 

r2 <- bmi ~ age

system <- list( r1, r2 )

imp <- mice(nhanes2, m = 5)
  m=5
  completed=lapply(1:5,function(i)complete(imp,i))
  fit.model <- systemfit(system, data= completed[[1]])

上面为每个估算数据集生成了完整的系统拟合输出。这很棒，但我剩下的任务是汇集 SUR 生成的整个输出。

我也未能成功尝试在 zelig 中运行该函数:

  completed.mi=do.call(Zelig:mi,completed)
  system=list(r1= bmi ~ hyp,r2=bmi~age)
  z.out=zelig(formula= system,model="sur",data=completed.mi)
  >Error: sur is not a supported model type.

最后，调用长形式的估算数据会产生很大的自由度，这并不反射(reflect)每个估算数据集中的实际案例数量(不包括示例)

我的问题是:

1) systemfit 软件包是否支持 MI 数据的 SUR？

2)如果是这样，您是否能够汇集所有估算数据集的输出？

3) 是否有替代包选项(除了 systemfit)用于在 R 中完成 SUR？

4) 如果 3 是否定的，是否有类似的分析可以实现相同的目标，并且是否有不同的包(例如 rms)可以支持 MI 数据池？

Answer 1

我认为小鼠不支持合并 SUR 的结果。您必须使用鲁宾规则手动汇总结果。我可以使用您的示例到达某个点，也许您可​​以从那里开始。

library(systemfit)
library(mice)
nhanes2

# add two imputation as example
imp <- mice(nhanes2, m = 2)
m=2

# create a data set with all the 
# complete imputed data sets
df<-mice::complete(imp, action="long", include=FALSE)

#create separate data frames for each mi
for(i in (df$.imp)) {
  nam <- paste0("df_", i)
  assign(nam, df[df$.imp==i,])
}

# Store the coefficients and se of each 
# sur in the M imputed data sets

M <-2 # number of imputations
M2 <- M*2 #number of total sur regressions
results <- as.data.frame(matrix(NA, nrow=M2, ncol = 4)) # will store here coefficients and se
colnames(results)<-c("coef_r1", "coef_r2", "se_r1", "se_r2")

# perform sur 
r1 <- bmi ~ hyp 
r2 <- bmi ~ age
system <- list( r1, r2 )

# start with first data set
fitsur1 <- systemfit(list( r1= r1, r2 = r2),
                data=df_1)
# start with second data set
fitsur2 <- systemfit(list( r1= r1, r2 = r2),
                 data=df_2)

# this can be warped in a loop
# but I could not do it...

# Extract coefficients
# Note I extract the coefficient 
# from only one age-group of r2, 
# Use same approach for the other
# extract coef from fitsur1
a <- as.data.frame(summary(fitsur1 )$coefficients[2,1])
colnames(a)<-c("coef_r1")
b <- as.data.frame(summary(fitsur1 )$coefficients[4,1])
colnames(b)<-c("coef_r2")
ab <- cbind(a,b)

# extract coef from fitsur2
c <- as.data.frame(summary(fitsur2 )$coefficients[2,1])
colnames(c)<-c("coef_r1")
d <- as.data.frame(summary(fitsur2 )$coefficients[4,1])
colnames(d)<-c("coef_r2")
cd <- cbind(c,d)


# Follow same approach to extract SE
# I cannot manage to extract them from 
# the 'fitsur' list ...

# merge extracted coef and se
results <- rbind(ab, cd)

# Then bind the coefficients and se
# from all imputed regressions

# Calculate the mean of pooled coefficients
pooled.coef_r1 <- mean(results$coef_r1)
pooled.coef_r2 <- mean(results$coef_r2)

计算池SE更加复杂 我用过这篇文章https://stats.stackexchange.com/questions/327237/calculating-pooled-p-values-manually

# example for se_r1
(betweenVar <- mean(results[,3])) # mean of variances
(withinVar <- sd(results[,1])^2) # variance of variances
(dfCorrection <- (nrow(results)+1)/(nrow(results))) # dfCorrection

(totVar <- betweenVar + withinVar*dfCorrection) # total variance
(pooledSE <- sqrt(totVar)) # standard error

我还没有研究过 p 值，但现在应该更容易