如何在R中使用predict函数来预测多年前拟合的逻辑回归模型？

Question

我有一个问题，一直在尝试解决但没有成功。已经搜索了两天以上，但仍然毫无头绪。如果答案已经存在而我没有找到，那么请见谅。

假设你有一个二元逻辑回归方程（二元模型）的回归方程，这是几年前估计出来的旧模型。因此，你知道参数βk（k = 1, 2，...，p），因为它们是过去估计出来的。但是你没有用于拟合模型的数据。

我的问题是：我是否可以将这个旧的估计逻辑模型作为一个对象引入R中（对应于逻辑回归模型）？

我想使用“预测”函数来证明这个逻辑回归模型与新数据集（现在的数据）相匹配，然后检查这个旧模型是否经得起时间的考验。而要使用这个函数，你需要逻辑回归模型的对象。

非常感谢您的帮助。

- perevales

---

这个问题似乎与统计有关，因此可能应该迁移到交叉验证。 - Thomas

用户正在尝试将他/她手头的一个方程式转换为一个对象。我认为这是一个相当适合在SO上提问的编程问题。 - Roman Luštrik

我可能会修改现有的模型，但那是作弊。 - Roman Luštrik

@RomanLuštrik 我猜，将方程转换为glm对象可能比基于底层统计数据从头构建算法更难。也许这仍然在范围内，我想。 - Thomas

Answer 1

根据我的评论，我认为你可以直接从系数计算预测结果。以下是一个示例，比较了predict.glm的输出和直接在数据上计算的预测概率：

# construct some data and model it
# y ~ x1 + x2
set.seed(1)
x1 <- runif(100)
x2 <- runif(100)
y <- rbinom(100,1,(x1+x2)/2)
data1 <- data.frame(x1=x1,x2=x2,y=y)
x3 <- runif(100)
x4 <- runif(100)
y2 <- rbinom(100,1,(x3+x4)/2)
data2 <- data.frame(x1=x3,x2=x4,y=y2)
glm1 <- glm(y~x1+x2,data=data1,family=binomial)

# extract coefs
#summary(glm1)
coef1 <- coef(glm1)

# calculate predicted probabilities for current data
tmp1 <- coef1[1] + (data1$x1*coef1[2]) + (data1$x2*coef1[3])
pr1 <- 1/(1+(1/exp(tmp1)))
# these match those from `predict`:
all.equal(pr1,predict(glm1,data1,type='response'))

# now apply to new data:
tmp2 <- coef1[1] + (data2$x1*coef1[2]) + (data2$x2*coef1[3])
pr2 <- 1/(1+(1/exp(tmp2)))
pr2

这显然不是一种通用的解决方案，也不能正确处理不确定性，但我认为这比篡改 predict 更好。

- Thomas

---

不错的设置 - 我会尝试运行 qqplot 将数据集与该模型生成的“模拟”数据进行比较。 - Carl Witthoft

Answer 2

您可以使用从已有系数创建的偏移量来创建glm拟合，然后使用常规的预测函数进行预测。例如，使用鸢尾花数据（首先在真实数据上拟合模型，然后使用虚拟数据和第一次拟合的系数来拟合新模型）：

fit1 <- glm( I(Species=='versicolor') ~ Petal.Length + Petal.Width, 
   data=iris, family=binomial )
coef(fit1)

dummydata <- data.frame( Petal.Length = rnorm(10), Petal.Width=rnorm(10),
    Species = rep(c('versicolor','other'), each=5) )

fit2 <- glm( I(Species=='versicolor') ~ 0 + 
  offset(-2.863708 + 1.563076*Petal.Length - 3.153165*Petal.Width),
    data=dummydata, family=binomial )

pred1 <- predict(fit1, newdata=iris)
pred2 <- predict(fit2, newdata=iris)
plot(pred1,pred2)
abline(0,1, col='green')

- Greg Snow