具有归一化数据池结果的R中看似不相关的回归

Question

我试图用R.中的systemfit软件包来完成看似不相关的回归(SUR)，然而，用多个推测的数据(用鼠标包)来完成这些分析并不简单。

在googling这个问题时，我发现有一个关于相同问题的被删除的帖子，它似乎使用了下面的例子(归功于海报，小编辑)。

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.

最后，调用计算数据的长形式会产生很大的自由度，它不反映每个估算数据集中的实际案例数(示例不包括在内)。

我的问题是：

1)系统是否支持MI数据？

2)如果是这样的话，您是否能够将输出集中在所有推测的数据集中？

( 3)在R中完成SUR是否有替代的一揽子方案(系统适配除外)？

4)如果3为否，是否有类似的分析将实现同样的目标，是否有不同的一揽子(如均方根)可能支持集中管理信息数据？

Answer 1

我不认为老鼠支持汇集SUR的结果。您必须使用Rubin的规则手动汇总结果。我可以用你的例子讲到一个特定的点，也许你可以从那里学到。

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更复杂--我使用了这个post 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值，但是现在应该更容易了。