首页 > 解决方案 > 在 R 中使用一组 beta 进行预测

问题描述

假设我有使用 MCMC 的 Beta 后验图。对于每一行,我都有一组测试版。是否可以将这一行 beta 转换为 R 中的模型对象,以便我可以使用 predict() 函数?具体来说,一些 beta 是针对分类随机变量的,所以如果我想手动应用 beta 会很困难。

我认为要手动完成,我必须将每个分类变量转换为多列指标变量。

标签: rregressionmcmc

解决方案


您可以通过构建模型矩阵并将其乘以系数向量来做到这一点。

假设这是我们的模型公式。Species是一个分类变量:

> m = Sepal.Length~Petal.Length+Petal.Width+Species

现在根据它和数据创建一个模型矩阵:

> mm = model.matrix(m, data=iris)
> head(mm)
  (Intercept) Petal.Length Petal.Width Speciesversicolor Speciesvirginica
1           1          1.4         0.2                 0                0
2           1          1.4         0.2                 0                0
3           1          1.3         0.2                 0                0
4           1          1.5         0.2                 0                0
5           1          1.4         0.2                 0                0
6           1          1.7         0.4                 0                0

具有Species三个级别,您可以看到前几行对于两个物种列都为零,因此它们属于其他物种。

要对一组 5 个系数进行预测,它们需要与行具有相同的顺序并具有相同的“对比”,即物种的编码使用两个虚拟变量以相同的方式完成。因此,您的“beta”值必须与您从lm数据中获得的系数看起来相同:

> mfit = lm(m, data=iris)
> mfit$coefficients
      (Intercept)      Petal.Length       Petal.Width Speciesversicolor 
      3.682982011       0.905945867      -0.005995401      -1.598361502 
 Speciesvirginica 
     -2.112646781 

现在您可以通过矩阵乘法获得任何一组系数的预测:

> head(mm %*% mfit$coefficients)
      [,1]
1 4.950107
2 4.950107
3 4.859513
4 5.040702
5 4.950107
6 5.220692

如果我做对了,这些值应该等于输出predict

> head(predict(mfit))
       1        2        3        4        5        6 
4.950107 4.950107 4.859513 5.040702 4.950107 5.220692 

如果你的系数不同,你应该通过矩阵乘法得到不同的预测:

> head(mm %*% c(0,1,.2,.2,.4))
  [,1]
1 1.44
2 1.44
3 1.34
4 1.54
5 1.44
6 1.78

推荐阅读