定义先验多项式回归。`MCMCglmm` 的案例研究

机器算法验证 r 贝叶斯 分类数据 马尔可夫链蒙特卡罗 事先的
2022-04-04 02:35:03

我不明白如何定义多项回归的先验分布。

  • 鉴于响应实际上没有单位(而只是类别),应该以什么单位设置先验概率?它应该是一个概率,还是一个赔率的对数?
  • 先验的维度是否取决于分类响应变量中的水平数?
  • 分类响应变量的方差 - 协方差矩阵的含义是什么?

这是我的数据的一个子集(变量名称和结果已被重命名)。处理这些数据的一个例子MCMCglmm会很棒!

df=read.table(text="y     x1       x2
    1  yellow 106.00 6.190476
    2  yellow 120.00 5.254762
    3  yellow  57.00 6.202381
    4  yellow 115.33 5.652381
    5  yellow 175.00 6.154762
    6  yellow  74.00 8.285714
    7  yellow 104.67 3.766667
    8  yellow  95.50 7.976190
    9  yellow 108.00 8.792857
    10 yellow 121.33 7.935714
    11 yellow  66.67 6.969048
    12 yellow  30.00 7.333333
    13 yellow  45.00 6.811905
    14 yellow  70.00 7.550000
    15 yellow  48.00 7.316667
    16 yellow 211.00 4.650000
    17 yellow  69.00 8.369048
    18 yellow 110.50 6.621429
    19 yellow 203.00 6.095238
    20 yellow  75.33 8.211905
    21 yellow 207.33 6.211905
    22 yellow  54.00 7.961905
    23 yellow  74.00 7.019048
    24 yellow 113.00 4.221429
    25 yellow  23.00 7.942857
    26 yellow  80.00 7.511905
    27 yellow 257.00 7.878571
    28 yellow 211.00 7.754762
    29 yellow  99.00 8.016667
    30 yellow 120.00 7.728571
    31 yellow 222.50 5.840476
    32 yellow  44.00 4.209524
    33 yellow  63.00 6.614286
    34 yellow  57.00 8.669048
    35 yellow 223.33 7.033333
    36 yellow 128.00 6.754762
    37 yellow 128.00 5.561905
    38 yellow 121.00 7.471429
    39 yellow  70.00 7.445238
    40 yellow  85.67 5.261905
    41 yellow 113.33 8.509524
    42 yellow  82.00 6.697619
    43    red 207.33 4.180952
    44    red 167.67 5.302381
    45    red 366.50 7.102381
    46    red 230.00 4.942857
    47    red 201.00 5.754762
    48    red 226.00 9.076190
    49    red 193.33 7.066667
    50    red 170.00 7.314286
    51    red 361.33 7.502381
    52   blue 154.00 4.342857
    53    red 199.33 6.361905
    54   blue  97.00 7.750000
    55   blue  82.33 6.209524
    56   blue  55.67 5.321429
    57   blue  47.50 5.911905
    58   blue  15.67 7.185714
    59   blue  96.50 6.452381
    60   blue 202.33 8.576190
    61   blue 157.00 6.669048
    62   blue 117.33 5.828571
    63   blue 105.67 8.485714
    64   blue 108.67 5.714286
    65   blue 296.67 5.852381
    66   blue 206.50 6.826190
    67   blue  88.50 6.178571
    68   blue 163.00 7.833333
    69   blue 151.50 8.983333")

这是一个MCMCglmm默认先验导致错误消息的调用

set.seed(12)
m = MCMCglmm(y ~  -1  + trait:(x1) + trait:(x2)  , rcov = ~ us(trait):units,
 data = df, family = "categorical", verbose = TRUE, burnin = 8000,
 nitt = 40000, thin = 50)

ill-conditioned G/R structure (CN = 24007848728601288.000000):
use proper priors if you haven't or rescale data if you have
1个回答

好问题!实际上,我已经尝试了相当长一段时间的想法,所以我将分享我的经验。我还是这个领域的新手,所以我希望我没有在符号上犯任何错误。

always,我的意思是always,将连续变量标准化为 mean=0 和 sd=1(甚至 sd=2)。查看 Andrew Gelman 的一些博客文章或文章。只是google Andrew Gelman 标准化,有很多好的论文和帖子。

将您的系数视为 log(odds-ratios) (参考参考类别,解释如下)。有关深入讨论,请参阅此答案Andrew Gelman 也有一些关于先验的建议,例如cauchynormal(0,1)他的论文是关于逻辑回归的,但我发现这些建议也扩展到多结果回归。

先验的维度确实取决于结果的数量。如果你有三个结果,你就有这三个线性依赖:

y1=β0,1+iβi,1xi,1y2=β0,2+iβi,2xi,2y3=β0,3+iβi,3xi,3

第二个下标表示结果。通常,您会对这些系数中的每一个进行先验。我将用截距来说明β0,k

β0,1N(0,1)β0,2N(0,1)β0,3N(0,1)

只是为了完成数学,结果的概率是,k

pk=ykkyk可能性是yMultinomial(p)

请注意,由于可识别性问题,您经常会强制将其中一个结果作为“参考”结果。Wiki有详细说明原因。实际上,这意味着您将参考类别的系数强制为 0,因此,因为preference=1exp(0)=1

我从来没有使用过 MCMCglmm,所以恐怕我无法回答任何具体的问题。