R中二項式glm響應的輸入格式
在
R
中,有三種方法可以使用glm
函數格式化邏輯回歸的輸入數據:
- 對於每個觀察,數據可以是“二進制”格式(例如,對於每個觀察,y = 0 或 1);
- 數據可以是“Wilkinson-Rogers”格式(例如,
y = cbind(success, failure)
),每一行代表一種治療;或者- 對於每個觀察,數據可以採用加權格式(例如,y = 0.3,權重 = 10)。
所有三種方法都產生相同的係數估計,但自由度和產生的偏差值和 AIC 分數不同。最後兩種方法具有較少的觀察值(因此具有較少的自由度),因為它們將每種處理用於觀察的數量,而第一種方法將每種觀察用於觀察的數量。
**我的問題:**使用一種輸入格式比另一種輸入格式有數字或統計優勢嗎?我看到的唯一優勢是不必重新格式化數據
R
以與模型一起使用。我查看了glm 文檔,在網絡上搜索,以及這個站點,發現了一個相關的帖子,但沒有關於這個主題的指導。
這是一個演示此行為的模擬示例:
# Write function to help simulate data drc4 <- function(x, b =1.0, c = 0, d = 1, e = 0){ (d - c)/ (1 + exp(-b * (log(x) - log(e)))) } # simulate long form of dataset nReps = 20 dfLong <- data.frame(dose = rep(seq(0, 10, by = 2), each = nReps)) dfLong$mortality <-rbinom(n = dim(dfLong)[1], size = 1, prob = drc4(dfLong$dose, b = 2, e = 5)) # aggregate to create short form of dataset dfShort <- aggregate(dfLong$mortality, by = list(dfLong$dose), FUN = sum) colnames(dfShort) <- c("dose", "mortality") dfShort$survival <- nReps - dfShort$mortality dfShort$nReps <- nReps dfShort$mortalityP <- dfShort$mortality / dfShort$nReps fitShort <- glm( cbind(mortality, survival) ~ dose, data = dfShort, family = "binomial") summary(fitShort) fitShortP <- glm( mortalityP ~ dose, data = dfShort, weights = nReps, family = "binomial") summary(fitShortP) fitLong <- glm( mortality ~ dose, data = dfLong, family = "binomial") summary(fitLong)
除了概念上的清晰性之外,沒有統計上的理由偏愛其中一個。儘管報告的偏差值不同,但這些差異完全是由於飽和模型造成的。因此,使用模型之間的相對偏差進行的任何比較都不會受到影響,因為飽和模型對數似然會抵消。
我認為通過顯式偏差計算很有用。
模型的偏差為 2*(LL(Saturated Model) - LL(Model))。假設你有 $ i $ 不同的細胞,其中 $ n_i $ 是單元格中的觀察數 $ i $ , $ p_i $ 是單元格中所有觀測值的模型預測 $ i $ , 和 $ y_{ij} $ 是觀測值(0 或 1) $ j $ - 細胞中的第一次觀察 $ i $ .
長表
(建議或空)模型的對數似然是 $$ \sum_i\sum_j\left(\log(p_i)y_{ij} + \log(1 - p_i)(1 - y_{ij})\right) $$
飽和模型的對數似然為$$ \sum_i\sum_j \left(\log(y_{ij})y_{ij} + \log(1 - y_{ij})(1-y_{ij})\right). $$這等於 0,因為 $ y_{ij} $ 為 0 或 1。 注意 $ \log(0) $ 未定義,但為方便起見,請閱讀 $ 0\log(0) $ 作為簡寫 $ \lim_{x \to 0^+}x\log(x) $ , 為 0。
短形式(加權)
請注意,二項式分佈實際上不能採用非整數值,但我們仍然可以通過使用每個單元格中觀察到的成功的分數作為響應來計算“對數似然”,並在對數似然計算中加權每個總和該單元格中的觀察次數。
$$ \sum_in_i \left(\log(p_i)\sum_jy_{ij}/n_i + \log(1 - p_i)(1 - \sum_j(y_{ij}/n_i)\right) $$
這完全等於我們上面計算的模型偏差,您可以通過將總和拉到 $ j $ 盡可能採用長式方程。
同時飽和偏差是不同的。由於我們不再有 0-1 響應,即使每次觀察只有一個參數,我們也不能準確地得到 0。相反,飽和模型對數似然是
$$ \sum_i n_i\left(\log(\sum_jy_{ij}/n_i)\sum_jy_{ij}/n_i + \log(1 - \sum_jy_{ij}/n_i)(1-\sum_jy_{ij}/n_i)\right). $$
在您的示例中,您可以驗證該數量的兩倍是兩個模型報告的空值和剩餘偏差值之間的差異。
ni = dfShort$nReps yavg = dfShort$mortalityP sum.terms <-ni*(log(yavg)*yavg + log(1 - yavg)*(1 - yavg)) # Need to handle NaN when yavg is exactly 0 sum.terms[1] <- log(1 - yavg[1])*(1 - yavg[1]) 2*sum(sum.terms) fitShortP$deviance - fitLong$deviance