我'm a bit lost on whether there'是一种更简单的方法，可以使用 predict 从更复杂的回归框架中得出预测结果 .

举个例子，如下：

NN<-1e4
data<-data.table(trt=sample(paste("Treatment",1:3),NN,T),
                 qtl=sample(paste0("Q",1:2),NN,T),
                 grp=sample(4,NN,T),
                 cat=sample(paste("Category",LETTERS[1:3]),NN,T),
                 val=rnorm(NN,10)^2)

data[,out:=140+5*(trt=="Treatment 2")+3*(trt=="Treatment 3")+
       8*(qtl=="Q2")-4*(trt=="Treatment 2"&qtl=="Q2")+
       7*(trt=="Treatment 3"&qtl=="Q2")-4*(grp==2)+
       6*(grp==3)-10*(grp==4)-6*(cat=="Category B")+
       2*(cat=="Category C")-1.8*val+rnorm(NN,10)>0]

llog<-glm(out~trt*qtl+as.factor(grp)+cat+val,data=data,family=binomial(link="logit"))

现在，我希望通过 trt ， qtl ， all other predictors held at sample averages ，以 out = 1的概率退出此预测 .

假设我想要这个，在这种情况下，在3 x 2矩阵（或表格等）中，行对应于 trt ，列为 qtl

特别是由于因子变量的存在，这很复杂 - “保持在样本平均值”意味着我们需要插入每组中观察的百分比，并且我不确定如何以干净的方式这样做 .

当然，真正漫长的方法是：

1）为“其他预测因子”设置一个样本平均值向量：

oth.avg<-c(1,unlist(data[,lapply(list(grp==2,grp==3,grp==4,cat=="Category B",
                                  cat=="Category C",val),mean)]))

2）乘以相应的系数

x.beta.oth<-sum(llog$coefficients[c(1,5:10)]*oth.avg)

3）为 trt ， qtl 和 trt x qtl 术语设置"matrices" :(我说"matrices"因为他们're conceptually drawn from underlying matrices, but it'更简洁，在一个维度中指定它们）

main.coef<-llog$coefficients[c(2:4,11:12)]
trt.mat<-rep(c(0,main.coef[1:2]),2)
qtl.mat<-rep(c(0,main.coef[3]),each=3)
tq.mat<-c(rep(0,3),rbind(rep(0,1),matrix(main.coef[4:5],ncol=1)))

（在这个小例子中，如此指定它们是过分的，但是当表格是4x4时，简约性开始显示出来 .

4）添加所有内容以获得预测的潜在指数

lat.pred<-trt.mat+qtl.mat+tq.mat+x.beta.oth

5）最后，通过逻辑公式将这些变换为预测概率：p /（1 p）

pred.prob<-matrix(exp(lat.pred)/(1+exp(lat.pred)),ncol=2,
              dimnames=list(paste("Treatment",1:3),c("Q1","Q2")))

                      Q1           Q2
Treatment 1 3.515452e-28 6.000161e-22
Treatment 2 1.989195e-24 1.519577e-27
Treatment 3 3.380796e-26 1.633097e-20

我错过了什么吗？我没有看到我可以输入 predict 以获得此输出...

r data.table glm predict