R 中看似无关的回归与估算的数据池结果

Question

我正在尝试使用 R 中的 systemfit 包完成看似无关的回归 (SUR)。但是，使用乘法插补数据（使用 mice 包）完成这些分析并不简单。

在谷歌搜索这个问题后，我看到关于相同问题的删除 post，它似乎使用了以下示例（归功于 poster，稍作修改）

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

以上为每个估算数据集生成完整的 systemfit 输出。这很好，但我的任务是汇集 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) systemfit 包是否支持 MI 数据的 SUR？

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

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

4) 如果 3 是否定的，是否有类似的分析可以实现相同的目标，是否有不同的包（例如，rms）可能支持 MI 数据的汇集？

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 值，但现在应该更容易了

R 中看似无关的回归与估算的数据池结果

Seemingly unrelated regression in R with imputed data-Pooling results

regression

r

pool

r-mice

imputation