從glm模型中提取參考類別的最佳方法？

我正在寫接受full和reducedglm對象概括爲感興趣varofint一個變量和相互作用可變interaction_var相互作用的結果（由full對象提取結果varofint有關執行lrtest並使用svycontrast的函數每個級別的interaction_var）。樣本數據：從glm模型中提取參考類別的最佳方法？

x <- data.frame(outcome=rbinom(100,1,.3),varofint=rnorm(100), interaction_var=sample(letters[1:3],100,replace=TRUE)) 

reduced <- glm(outcome~varofint+interaction_var,data=x) 
full <- glm(outcome~varofint*interaction_var,data=x)

我想知道提取參考類別說（full）GLM模型的最佳途徑。我能明顯這樣做

levels(full$data$interaction_var)[1]

，但是這會是一個「安全」的方法來提取給定輸入到contrasts參數參考類別？看起來，如果選擇選擇SAS對比度，則此方法可產生interactionv_var的等級，該等級不是在模型中用作參考類別的等級。以下會更安全嗎？

mf <- model.frame(full) 
setdiff(rownames(contrasts(mf[, "interaction_var"])), colnames(contrasts(mf[, "interaction_var"])))

或類似

names(which(apply(contrasts(mf[, "interaction_var"]),1,function(.v){all(.v==0)})))

我缺少一個更簡單的方法來提取參考類別？

來源

2013-12-23 Michael

如果沒有參考類別？參考類別僅適用於治療對比。 –

好吧，這是一個好的開始。所以，如果相應的交互變量沒有處理對比，那麼函數應該返回一個錯誤... – Michael

下面是這個任務的函數：

refCat <- function(model, var) { 
    cs <- attr(model.matrix(model), "contrasts")[[var]] 
    if (is.character(cs)) { 
    if (cs == "contr.treatment") 
     ref <- 1 
    else stop("No treatment contrast") 
    } 
    else { 
    zeroes <- !cs 
    ones <- cs == 1 
    stopifnot(all(zeroes | ones)) 
    cos <- colSums(ones) 
    stopifnot(all(cos == 1)) 
    ros <- rowSums(ones) 
    stopifnot(sum(!ros) == 1 && sum(ros) != ncol(cs)) 
    ref <- which(!ros) 
    } 
    return(levels(model$data[[var]])[ref]) 
}

的功能，如果變量var並不表示爲處理停止的對比。

例子：

refCat(reduced, "interaction_var") 
# [1] "a" 
refCat(full, "interaction_var") 
# [1] "a"

來源

2013-12-24 09:49:26

謝謝，這似乎工作。在「is.character（cs）」爲「FALSE」的情況下？ – Michael

@Michael'cs'可以是具有對比度函數名稱的字符串，例如'contr.treatment'，或數字對比矩陣。 –

從glm模型中提取參考類別的最佳方法？

回答

相關問題