自動比較鼠標glm.mids中的嵌套模型

Question

我有一個來自R的mice包的乘法 - 估算模型，其中有很多因子變量。例如：

library(mice) 
library(Hmisc) 

# turn all the variables into factors 
fake = nhanes 
fake$age = as.factor(nhanes$age) 
fake$bmi = cut2(nhanes$bmi, g=3) 
fake$chl = cut2(nhanes$chl, g=3) 

head(fake) 
    age   bmi hyp  chl 
1 1  <NA> NA  <NA> 
2 2 [20.4,25.5) 1 [187,206) 
3 1  <NA> 1 [187,206) 
4 3  <NA> NA  <NA> 
5 1 [20.4,25.5) 1 [113,187) 
6 3  <NA> NA [113,187) 

imput = mice(nhanes) 

# big model 
fit1 = glm.mids((hyp==2) ~ age + bmi + chl, data=imput, family = binomial) 

我想通過針對在滴一個變量中的每個可能的嵌套模型中測試的完整模型，以測試在所述模型（未每個級別的指示器變量）每個整個因子變量的意義一次。手動，我可以這樣做：

# small model (no chl) 
fit2 = glm.mids((hyp==2) ~ age + bmi, data=imput, family = binomial) 

# extract p-value from pool.compare 
pool.compare(fit1, fit2)$pvalue 

如何自動爲我的模型中的所有因子做這些事情？建議drop1非常有幫助的功能a previous question - 現在我想要做的事情完全一樣，除了mice情況。

可能有用注：的pool.compare一個惱人的特點是，它似乎要在更大的模型「額外」的變量被放置與該小模型共享的那些之後。

Answer 1

在按照pool.compare所需的順序排列它們之後，可以使用循環遍歷預測變量的不同組合。

因此，使用從上面的數據fake - 調整了類別

library(mice) 
library(Hmisc) 
# turn all the variables into factors 
# turn all the variables into factors 
fake <- nhanes 
fake$age <- as.factor(nhanes$age) 
fake$bmi <- cut2(nhanes$bmi, g=2) 
fake$chl <- cut2(nhanes$chl, g=2) 

# Impute 
imput <- mice(fake, seed=1) 

# Create models 
# - reduced models with one variable removed 
# - full models with extra variables at end of expression 
vars <- c("age", "bmi", "chl") 

red <- combn(vars, length(vars)-1 , simplify=FALSE) 
diffs <- lapply(red, function(i) setdiff(vars, i)) 
(full <- lapply(1:length(red), function(i) 
          paste(c(red[[i]], diffs[[i]]), collapse=" + "))) 
#[[1]] 
#[1] "age + bmi + chl" 

#[[2]] 
#[1] "age + chl + bmi" 

#[[3]] 
#[1] "bmi + chl + age" 

(red <- combn(vars, length(vars)-1 , FUN=paste, collapse=" + ")) 
#[1] "age + bmi" "age + chl" "bmi + chl" 

該機型現在在正確的順序傳遞給glm來電的號碼。我也換成glm.mids方法，因爲它已取代with.mids - 見?glm.mids

out <- vector("list", length(red)) 

for(i in 1:length(red)) { 

    redMod <- with(imput, 
       glm(formula(paste("(hyp==2) ~ ", red[[i]])), family = binomial)) 

    fullMod <- with(imput, 
       glm(formula(paste("(hyp==2) ~ ", full[[i]])), family = binomial)) 

    out[[i]] <- list(predictors = diffs[[i]], 
        pval = c(pool.compare(fullMod, redMod)$pvalue)) 
    } 

do.call(rbind.data.frame, out) 
# predictors  pval 
#2   chl 0.9976629 
#21  bmi 0.9985028 
#3   age 0.9815831 

# Check manually by leaving out chl 
mod1 <- with(imput, glm((hyp==2) ~ age + bmi + chl , family = binomial)) 
mod2 <- with(imput, glm((hyp==2) ~ age + bmi , family = binomial)) 
pool.compare(mod1, mod2)$pvalue 
#   [,1] 
#[1,] 0.9976629 

你會使用這個數據集得到了很多警告

編輯

你可以在一個函數包裝這個

impGlmDrop1 <- function(vars, outcome, Data=imput, Family="binomial") 
{ 

    red <- combn(vars, length(vars)-1 , simplify=FALSE) 
    diffs <- lapply(red, function(i) setdiff(vars, i)) 
    full <- lapply(1:length(red), function(i) 
         paste(c(red[[i]], diffs[[i]]), collapse=" + ")) 
    red <- combn(vars, length(vars)-1 , FUN=paste, collapse=" + ") 

    out <- vector("list", length(red)) 
    for(i in 1:length(red)) { 

    redMod <- with(Data, 
       glm(formula(paste(outcome, red[[i]], sep="~")), family = Family)) 
    fullMod <- with(Data, 
       glm(formula(paste(outcome, full[[i]], sep="~")), family = Family)) 
    out[[i]] <- list(predictors = diffs[[i]], 
        pval = c(pool.compare(fullMod, redMod)$pvalue) ) 
    } 
    do.call(rbind.data.frame, out) 
} 

# Run 
impGlmDrop1(c("age", "bmi", "chl"), "(hyp==2)")