使用 glm.fit 后的方差分析

Question

我想执行似然比检验以确定 DOE 中模型项的功效。到目前为止，我一直在使用 glm fit 的 p 值来执行此操作，并且一切都很好。当我开始使用 anova 函数时，我意识到似乎没有设计为接受 glm.fit 函数输入的 anova 函数，只有 glm 函数。这是我想做的一个例子：

X # This is a model matrix from matrix.model
y # These are the y values for the fit
tfit = glm.fit(X, y, family = poisson())
anova(tfit, test = 'LRT')

通常我会假设anova 函数调用只需更改为anova.glm，但事实并非如此。如何使glm.fit 函数输出与anova 函数输入兼容？

Answer 1

问题在于glm.fit 输出的不是glm 类，而是包含有关模型的各种数据的原始列表。这不能提供给anova.glm，因为这个函数需要glm 函数产生的glm 类的对象。如果您有可用的原始数据（因此未转化为模型矩阵，则可以将glm 函数应用于此以产生所需的结果。

X <- matrix(c(runif(10), rnorm(10)), ncol = 2)
y <- round(runif(10, 1, 5))
X.mm <- model.matrix(y ~ X)
model.fit.1 <- glm.fit(X.mm, y, family = poisson())
class(model.fit.1)
model.fit.2 <- glm(y ~ X, family = "poisson")
class(model.fit.2)

anova(model.fit.2, test = "LRT")

如果您不能使用glm 函数而必须使用glm.fit，那么您可以从glm.fit 输出中自己构建LRT。首先使用以下函数

LRT.glm.fit <- function(glm.fit.mod){
  df.null <- glm.fit.mod$df.null
  df.mod <- glm.fit.mod$df.residual
  dev.null <- glm.fit.mod$null.deviance
  dev.mod <- glm.fit.mod$deviance
  dev.diff <- dev.null - dev.mod
  p.value <- 1 - pchisq(dev.null - dev.mod, df.null - df.mod)
  output <- c(round(df.null), round(df.mod), dev.null, dev.mod, p.value)
  names(output) <- c("df.null", "df.mod", "dev.null", "dev.mod", "p.value")
  output
}