使用RJags的Dirichlet先验分类数据模型

xxhby3vn  于 2023-04-18  发布在  其他
关注(0)|答案(1)|浏览(72)

我试图用jags写一个模型。问题是我对jags语言真的很陌生。我有一个数据集，基本上是一个调查，有6个问题和一个预测变量。所有都是分类的，其中一些有3个可能的答案，其他有2个，（预测有3个可能的答案）。为了使用锯齿创建贝叶斯模型，我决定使用分类分布来为Y（可能性）和Dirichlet（先验）建模。
因为我在X和Y中都有缺失值。在编写了一个模型来包含和计算这些值之后，R给了我一个致命的错误，我不能运行模型。我在网上找不到任何这样的案例/模型的来源来尝试和复制模型到我的案例中。无论如何，下面我有我的模型，所以请如果你发现有什么奇怪的，请告诉我，所以我修复的问题，如果一切都与模型似乎很好，你知道是什么原因导致的错误，我也很高兴听到这个。
产品规格：- 操作系统：MacOS Intel Core i5 - R StudioR 4.1.2 - rjags版本：4-13

Model : 
model {
  # Likelihood
 for (i in 1:209) {
  FIN[i] ~ dcat(p[i,1:3])
  logit(p[i,1]) <- alpha[1] + beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
  logit(p[i,2]) <- alpha[2] + beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
  logit(p[i,3]) <- alpha[3] + beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
 
  # -- Add model for dependent variables --
  AGE[i] ~ dcat(p_age)
  SEX[i] ~ dbern(p_sex)
  INC[i] ~ dcat(p_inc)
  POL[i] ~ dcat(p_pol)
  STAT[i]~ dbern(p_stat)
 }
  # Priors
  for (j in 1:6) {
  beta[j] ~ dnorm(0, 0.001) # Normal prior for beta
 }
  for (k in 1:3) {
  alpha[k] ~ dnorm(0, 0.001) # Normal prior for alpha
  }
  # -- Priors for dependent variables --
  p_age ~ ddirch(c(0.3444976,0.3157895,0.3349282))
  p_inc ~ ddirch(c(0.2248804,0.3971292,0.2870813))
  p_pol ~ ddirch(c(0.4019139,0.1913876,0.3732057 ))
  p_sex ~ dbeta(2,1)
  p_stat ~ dbeta(2,1)

如果模型中的任何修改可能适用于我的数据情况，请随意提出。
我的前10行数据如下：

AGE SEX INC POL AREA FIN STAT
1    3   0   2   2    1   2    1
2    2   0   3   3    1   3    1
3    1   0   2  NA    1   2    1
4    3   1   2   1    1   1    0
5    3   1   3   3    3   2   NA
6    1   0   2   1    3   3   NA
7    3   1   1   3    3   1    1
8    1   0   1   3    2   1    0
9    3   1  NA   3    3   2    0
10   1   0  NA   1    1   2    1

预测变量为FIN，其余为响应变量。

r

1条答案

按热度按时间
2uluyalo

2uluyalo1#

看起来你的似然是多项式和有序似然的混合。在多项式似然中，你会为结果类别中的所有变量设置不同的系数，但参考组的所有系数都将设置为零。对于有序logit似然，你会像你在这里做的那样-所有系数都被限制为在结果类别中相同，除了截距。然而，截距被强制排序，并且从排序的截点计算的概率是累积的（而不是针对每个类别）。因为我不确定你的意图，让我们看看两者。首先，多项式似然：

for (i in 1:209) {
    FIN[i] ~ dcat(p[i,1:3])
    log(q[i,1]) <- alpha[1] + beta[1, 1]*AGE[i] + beta[2, 1]*SEX[i] + beta[3, 1]*INC[i] + beta[4, 1]*POL[i] + beta[5, 1]*AREA[i] + beta[6, 1]*STAT[i]
    log(q[i,2]) <- alpha[2] + beta[1, 2]*AGE[i] + beta[2, 2]*SEX[i] + beta[3, 2]*INC[i] + beta[4, 2]*POL[i] + beta[5, 2]*AREA[i] + beta[6, 2]*STAT[i]
    log(q[i,3]) <- alpha[3] + beta[1, 3]*AGE[i] + beta[2, 3]*SEX[i] + beta[3, 3]*INC[i] + beta[4, 3]*POL[i] + beta[5, 3]*AREA[i] + beta[6, 3]*STAT[i]
    for(k in 1:3){
      p[i,k] <- q[i,k]/sum(q[i,])
    }
  }
  for (j in 1:6) {
    beta[k,1] <- 0
    for(k in 2:3){
      beta[j,k] ~ dnorm(0, 0.001) # Normal prior for beta
    }
  }
  alpha[1] <- 0
  for (k in 2:3) {
    alpha[k] ~ dnorm(0, 0.001) # Normal prior for alpha
  }

注意，我们将logit(p[i,k])改为log(q[i,k])，因为观察i在类别k中的概率为

我们将第一类系数设置为零作为参考。
有序模型似然性如下所示：

for (i in 1:209) {
  FIN[i] ~ dcat(p[i,1:3])
  logit(q[i,1]) <- alpha[1] - beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
  logit(q[i,2]) <- alpha[2] - beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
  p[i,1] <- q[i,1]
  p[i,2] <- q[i,2]-q[i,1]
  p[i,3] <- 1-q[i,2]
}
for (j in 1:6) {
    beta[j] ~ dnorm(0, 0.001) # Normal prior for beta
}
for (k in 1:2) {
  a[k] ~ dnorm(0, 0.001) # Normal prior for alpha
}
alpha <- sort(a)

这里，alpha参数是截断点-假设alpha 0 = -Inf且alpha 3 = Inf。重要的一点是alpha参数的大小必须增加。
让我们尝试运行两个模型：

多项式模型

mnl_mod <- "model {
  # Likelihood
  for (i in 1:10) {
    FIN[i] ~ dcat(p[i,1:3])
    log(q[i,1]) <- alpha[1] + beta[1, 1]*AGE[i] + beta[2, 1]*SEX[i] + beta[3, 1]*INC[i] + beta[4, 1]*POL[i] + beta[5, 1]*AREA[i] + beta[6, 1]*STAT[i]
    log(q[i,2]) <- alpha[2] + beta[1, 2]*AGE[i] + beta[2, 2]*SEX[i] + beta[3, 2]*INC[i] + beta[4, 2]*POL[i] + beta[5, 2]*AREA[i] + beta[6, 2]*STAT[i]
    log(q[i,3]) <- alpha[3] + beta[1, 3]*AGE[i] + beta[2, 3]*SEX[i] + beta[3, 3]*INC[i] + beta[4, 3]*POL[i] + beta[5, 3]*AREA[i] + beta[6, 3]*STAT[i]
    for(k in 1:3){
      p[i,k] <- q[i,k]/sum(q[i,])
    }
        # -- Add model for dependent variables --
    AGE[i] ~ dcat(p_age)
    SEX[i] ~ dbern(p_sex)
    INC[i] ~ dcat(p_inc)
    POL[i] ~ dcat(p_pol)
    STAT[i]~ dbern(p_stat)
  }
  for (j in 1:6) {
    beta[j,1] <- 0
    for(k in 2:3){
      beta[j,k] ~ dnorm(0, .1) # Normal prior for beta
    }
  }
  alpha[1] <- 0
  for (k in 2:3) {
    alpha[k] ~ dnorm(0, .1) # Normal prior for alpha
  }
  # Priors
  # -- Priors for dependent variables --
  p_age ~ ddirch(c(0.3444976,0.3157895,0.3349282))
  p_inc ~ ddirch(c(0.2248804,0.3971292,0.2870813))
  p_pol ~ ddirch(c(0.4019139,0.1913876,0.3732057 ))
  p_sex ~ dbeta(2,1)
  p_stat ~ dbeta(2,1)
}"
dat <- read.table(text="AGE SEX INC POL AREA FIN STAT
1    3   0   2   2    1   2    1
2    2   0   3   3    1   3    1
3    1   0   2  NA    1   2    1
4    3   1   2   1    1   1    0
5    3   1   3   3    3   2   NA
6    1   0   2   1    3   3   NA
7    3   1   1   3    3   1    1
8    1   0   1   3    2   1    0
9    3   1  NA   3    3   2    0
10   1   0  NA   1    1   2    1", header=TRUE)

dl <- as.list(dat)
library(rjags)
#> Loading required package: coda
#> Linked to JAGS 4.3.1
#> Loaded modules: basemod,bugs
mnl_j <- jags.model(textConnection(mnl_mod), data=dl, n.chains=2)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 55
#>    Unobserved stochastic nodes: 24
#>    Total graph size: 268
#> 
#> Initializing model
update(mnl_j, 10000)
mnl_samps <- coda.samples(mnl_j, variable.names=c("alpha", "beta"), n.iter=5000)

型号总结：

summary(mnl_samps)
#> 
#> Iterations = 11001:16000
#> Thinning interval = 1 
#> Number of chains = 2 
#> Sample size per chain = 5000 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>               Mean    SD Naive SE Time-series SE
#> alpha[1]   0.00000 0.000  0.00000        0.00000
#> alpha[2]  -0.98913 2.607  0.02607        0.12603
#> alpha[3]  -1.44130 2.869  0.02869        0.12171
#> beta[1,1]  0.00000 0.000  0.00000        0.00000
#> beta[2,1]  0.00000 0.000  0.00000        0.00000
#> beta[3,1]  0.00000 0.000  0.00000        0.00000
#> beta[4,1]  0.00000 0.000  0.00000        0.00000
#> beta[5,1]  0.00000 0.000  0.00000        0.00000
#> beta[6,1]  0.00000 0.000  0.00000        0.00000
#> beta[1,2] -0.64053 1.351  0.01351        0.07672
#> beta[2,2] -0.45833 2.317  0.02317        0.07423
#> beta[3,2]  2.15004 1.633  0.01633        0.10340
#> beta[4,2] -0.52331 1.455  0.01455        0.09729
#> beta[5,2] -0.17880 1.503  0.01503        0.08903
#> beta[6,2]  1.86677 2.012  0.02012        0.06045
#> beta[1,3] -2.50116 1.827  0.01827        0.08436
#> beta[2,3] -2.66265 2.719  0.02719        0.06562
#> beta[3,3]  3.98474 2.097  0.02097        0.12879
#> beta[4,3] -0.86416 1.609  0.01609        0.08633
#> beta[5,3]  0.07456 1.597  0.01597        0.07076
#> beta[6,3]  0.45412 2.697  0.02697        0.09375
#> 
#> 2. Quantiles for each variable:
#> 
#>              2.5%     25%     50%     75%  97.5%
#> alpha[1]   0.0000  0.0000  0.0000  0.0000 0.0000
#> alpha[2]  -6.0760 -2.7514 -0.9464  0.7636 4.0737
#> alpha[3]  -7.2179 -3.3325 -1.3085  0.4784 3.9364
#> beta[1,1]  0.0000  0.0000  0.0000  0.0000 0.0000
#> beta[2,1]  0.0000  0.0000  0.0000  0.0000 0.0000
#> beta[3,1]  0.0000  0.0000  0.0000  0.0000 0.0000
#> beta[4,1]  0.0000  0.0000  0.0000  0.0000 0.0000
#> beta[5,1]  0.0000  0.0000  0.0000  0.0000 0.0000
#> beta[6,1]  0.0000  0.0000  0.0000  0.0000 0.0000
#> beta[1,2] -3.2657 -1.5539 -0.6686  0.2970 1.9910
#> beta[2,2] -5.0281 -2.0270 -0.4329  1.1616 4.0349
#> beta[3,2] -0.9553  1.0735  2.0957  3.1838 5.5768
#> beta[4,2] -3.5925 -1.4003 -0.4761  0.4134 2.3058
#> beta[5,2] -3.2009 -1.1905 -0.1978  0.8687 2.7151
#> beta[6,2] -2.1637  0.5484  1.8749  3.1853 5.8102
#> beta[1,3] -6.4089 -3.6391 -2.3866 -1.2600 0.8229
#> beta[2,3] -7.9316 -4.4922 -2.6946 -0.8260 2.7573
#> beta[3,3]  0.1021  2.5416  3.8960  5.3833 8.2670
#> beta[4,3] -4.0675 -1.9166 -0.8612  0.2320 2.2527
#> beta[5,3] -3.0094 -1.0119  0.0389  1.1190 3.3034
#> beta[6,3] -4.8934 -1.3242  0.4582  2.2557 5.8092

诊断：

novar <- which(apply(mnl_samps[[1]], 2, sd) == 0)
gelman.diag(mnl_samps[,-novar])
#> Potential scale reduction factors:
#> 
#>           Point est. Upper C.I.
#> alpha[2]        1.00       1.01
#> alpha[3]        1.00       1.00
#> beta[1,2]       1.01       1.04
#> beta[2,2]       1.00       1.02
#> beta[3,2]       1.00       1.01
#> beta[4,2]       1.00       1.01
#> beta[5,2]       1.00       1.01
#> beta[6,2]       1.00       1.00
#> beta[1,3]       1.00       1.02
#> beta[2,3]       1.00       1.00
#> beta[3,3]       1.00       1.01
#> beta[4,3]       1.01       1.01
#> beta[5,3]       1.00       1.01
#> beta[6,3]       1.00       1.01
#> 
#> Multivariate psrf
#> 
#> 1.01

有序Logit模型

ol_mod <- "model {
  # Likelihood
  for (i in 1:10) {
    FIN[i] ~ dcat(p[i,1:3])
    logit(q[i,1]) <- alpha[1] - beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
    logit(q[i,2]) <- alpha[2] - beta[1]*AGE[i] + beta[2]*SEX[i] + beta[3]*INC[i] + beta[4]*POL[i] + beta[5]*AREA[i] + beta[6]*STAT[i]
    p[i,1] <- q[i,1]
    p[i,2] <- q[i,2]-q[i,1]
    p[i,3] <- 1-q[i,2]
    AGE[i] ~ dcat(p_age)
    SEX[i] ~ dbern(p_sex)
    INC[i] ~ dcat(p_inc)
    POL[i] ~ dcat(p_pol)
    STAT[i]~ dbern(p_stat)
  }
  for (j in 1:6) {
      beta[j] ~ dnorm(0, 0.25) # Normal prior for beta
  }
  a[1] ~ dnorm(-1, 0.25) # Normal prior for alpha
  a[2] ~ dnorm(1, 0.25) # Normal prior for alpha
  alpha <- sort(a)
        # -- Add model for dependent variables --
  # Priors
  # -- Priors for dependent variables --
  p_age ~ ddirch(c(0.3444976,0.3157895,0.3349282))
  p_inc ~ ddirch(c(0.2248804,0.3971292,0.2870813))
  p_pol ~ ddirch(c(0.4019139,0.1913876,0.3732057 ))
  p_sex ~ dbeta(2,1)
  p_stat ~ dbeta(2,1)
}"
ol_j <- jags.model(textConnection(ol_mod), data=dl, n.chains=2, inits=list(a = c(-1,1)))
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 55
#>    Unobserved stochastic nodes: 18
#>    Total graph size: 201
#> 
#> Initializing model
update(ol_j, 10000)
ol_samps <- coda.samples(ol_j, variable.names=c("alpha", "beta"), n.iter=5000)

型号总结：

summary(ol_samps)
#> 
#> Iterations = 11001:16000
#> Thinning interval = 1 
#> Number of chains = 2 
#> Sample size per chain = 5000 
#> 
#> 1. Empirical mean and standard deviation for each variable,
#>    plus standard error of the mean:
#> 
#>             Mean     SD Naive SE Time-series SE
#> alpha[1] -1.0566 1.4416 0.014416        0.04071
#> alpha[2]  2.9102 1.4574 0.014574        0.04246
#> beta[1]  -1.1004 0.9484 0.009484        0.04393
#> beta[2]   1.4218 1.5966 0.015966        0.04189
#> beta[3]  -2.1456 1.1222 0.011222        0.05998
#> beta[4]   0.8146 0.8649 0.008649        0.03884
#> beta[5]  -0.3696 0.8421 0.008421        0.03201
#> beta[6]  -0.5356 1.3743 0.013743        0.03513
#> 
#> 2. Quantiles for each variable:
#> 
#>              2.5%     25%     50%      75%   97.5%
#> alpha[1] -3.91901 -2.0171 -1.0359 -0.06762  1.6986
#> alpha[2]  0.06765  1.9184  2.9027  3.88113  5.8073
#> beta[1]  -3.09245 -1.7250 -1.0667 -0.44667  0.6355
#> beta[2]  -1.69368  0.3684  1.4133  2.49551  4.5652
#> beta[3]  -4.50615 -2.8821 -2.0831 -1.35712 -0.1230
#> beta[4]  -0.88252  0.2309  0.8094  1.38204  2.5435
#> beta[5]  -2.08284 -0.9168 -0.3691  0.20091  1.2354
#> beta[6]  -3.19611 -1.4704 -0.5316  0.35694  2.1961

诊断：

gelman.diag(ol_samps)
#> Potential scale reduction factors:
#> 
#>          Point est. Upper C.I.
#> alpha[1]       1.00       1.01
#> alpha[2]       1.00       1.01
#> beta[1]        1.01       1.01
#> beta[2]        1.00       1.02
#> beta[3]        1.01       1.01
#> beta[4]        1.01       1.06
#> beta[5]        1.01       1.05
#> beta[6]        1.00       1.00
#> 
#> Multivariate psrf
#> 
#> 1.02

创建于2023-04-09使用reprex v2.0.2
注意，在上面的代码中，我增加了系数的精度，以获得看起来接近收敛的结果，但如果你有很多数据，你应该能够在你的真实的模型中放松这些数据。

赞(0）分享  回复(0）举报  2023-04-18
我来回答

相关问题

    查看更多

    热门标签

    更多
    JavaquerypythonNode开发语言requestUtil数据库Table后端算法LoggerMessageElementParser

    最新问答

    更多
    • 蜀ICP备13028337号-1 大数据知识库 https://www.saoniuhuo.com © All rights reserved
    • 本站内容来源互联网,如果侵犯您的权益请联系我们删除, 联系方式：448109455@qq.com