在R中使用任意系数的predict()函数

Question

我有一些由非R用户设置的逻辑回归模型系数。我想将这些系数导入R并在同一数据集上生成一些拟合度估计（ROC和混淆矩阵），与我的模型进行对比。我的第一个想法是使用以下方法将系数强制转换为现有的GLM对象： summary(fit)$coefficients[,1] <- y 或者 summary(fit)$coefficients <- x 其中y和x是包含我试图用于预测的系数的矩阵，而fit是预先创建的虚拟glm对象，用于数据集的拟合。当然，这只给我带来了错误。

是否有任何方法可以将任意系数向量传递给predict()函数或在模型中指定系数？我是否可以通过将一个向量传递到GLM的offset参数中来强制执行此操作？谢谢

编辑：如评论中所述，使用任意系数没有太多的统计基础。我有一个商业伙伴，他/她相信自己“知道”正确的系数，我正在尝试量化基于这些估计值与由正确模型生成的系数之间的预测能力损失。

编辑2：根据BondedDust的答案，我能够强制转换系数，但无法清除由于强制转换而返回的predict()的错误消息，似乎调用predict的predict.lm也查看系数的秩，这导致了错误。

- Stencil

---

作为对这个问题的回应，我创建了一个makeglm()函数，听起来在这种情况下会很有用。如果您提供一个更具体的样本，我们可能可以尝试一下。 - MrFlick

@Stencill 你可以通过手动将数据中的相关变量乘以给定的系数来计算预测值吗？例如，coefVector %*% t(cbind(1, dataVariables))。（其中coefVector是提供的系数向量，datavariables是与相关系数对应的数据） - user20650

@BondedDust 抱歉，上周末离开电脑了--从现在开始我会及时更新的。 - Stencil

@user20650 这似乎是最简单的解决方案。但是，我需要将一些因素转换为虚拟变量。 - Stencil

没错 - model.matrix 可以使这个过程变得简单直接。 - user20650

Answer 1

如果您跟踪通过predict.glm传递对象到predict.lm的代码，似乎需要更改的模型列表节点确实是fit$coefficients。然而，更改summary()对象将没有任何效果。在glm和lm对象中的[['coefficients']]不是像summary所生成的具有“Estimate”，“Std. Error”，“t value”，“Pr（> | t |）”列的矩阵，而只是系数的向量。

 fit$coefficients <- y
 newpred <- predict(fit)

如果您需要进一步使用 fit，可以复制并进行操作。

- IRTFM

---

这个可行。非常感谢。这有点像是一个hack，因为glm对象的其余部分（Pvalues等）现在不匹配了，但我成功地通过predict()传递了它。 - Stencil

取消那个，似乎我在将fit传递到predict中时出现了错误。>pred <- predict(fit, newdata = sample1)其中fit是包含修改系数的glm对象，返回以下错误：

Error in [.data.frame(beta, piv) : undefined columns selected 此外：警告信息： 在predict.lm（object，newdata，se.fit，scale = 1，type = ifelse（type == ：prediction from a rank-deficient fit may be misleading

- Stencil

这绝对是一个hack，但这绝对是你要求的。你使用的方法没有真正的统计依据。我不知道你在解决这个问题中的进展。你应该 A) 首先发布一个数据示例，B) 使用编辑过程更新你的问题，而不是在我的答案下发表评论。 - IRTFM

Answer 2

这不是对您发布的问题的答案 - BondedDust已经回答了 - 但是描述了一种替代方法，可以帮助您自己计算预测概率。

# Use the mtcars dataset for a minimum worked example
data(mtcars)

# Run a logistic regression and get predictions 
mod <- glm(vs ~ mpg + factor(gear) + factor(am), mtcars, family="binomial")
p1 <- predict(mod, type="response")

# Calculate predicted probabilities manually
m <- model.matrix(~ mpg + factor(gear) + factor(am), mtcars)[,]
p2 <- coef(mod) %*% t(m)
p2 <- plogis(p2)

all(p1 == p2)
#identical(as.numeric(p1), as.numeric(p2))

你可以用给你的系数向量替换 coef(mod)。model.matrix 会生成计算所需的虚拟变量 - 请检查顺序与系数向量相同。

- user20650

---

如果我想基于 probit 模型进行预测（即在 glm 命令中使用 family = binomial(link = "probit")），我该如何更改上面的内容（即 plogis(p2)）？ - rp1

@rp1；我不确定 - probit链接的公式是什么？（因为 plogis == 1/(1+e(-xb))）。我认为您将使用 pnorm 函数(pnorm(p2))，但最好在 https://stats.stackexchange.com/questions 上询问。 - user20650

@rp1; 好的，我刚刚运行了示例，将链接更改为 probit: family=binomial("probit")，然后 p2 <- pnorm(p2)，接着 all.equal(p1, as.numeric(p2), check.attributes = FALSE)，所以似乎 pnorm 是正确的方法。 - user20650

@user20650：非常感谢。我也进行了一些测试，可以确认p2 <- pnorm(p2)是正确的方法。感谢您的快速回复！ - rp1

Answer 3

或者，您可以使用类似以下的方法：

fit <- lm(Y ~ A + B + C, data=fakedata)

fit$coefficients <- c(1, 2, 3) # 这会将A、B、C的系数分别更改为1、2和3。

Y_hat_new <- predict(fit, new_fakedata) # 基于新的系数和/或新数据，计算出预测结果Y_hat_new。

如果按照模型矩阵的方法进行，结果应该是相同的。

- Runfirst