R语言jags例子

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

mf2 <- model.frame(formula_2a,study)
mt2 <- attr(mf2, "terms", exact = TRUE)
x2_a <- as.data.frame(model.matrix(mt2, mf2))[,-1]
y2_a <- model.response(mf2)
n2_a = length(y2_a)
model_code_2a = function(){
  # Likelihood
  for (i in 1:n) {
    y[i] ~ dnorm(alpha + genderOtra * x1[i] + genderWoman * x2[i] + on_a_dietYes * x3[i] +
                         alcoholYes * x4[i] + drugsYes * x5[i] + sickYes * x6[i] +
                         other_factorsYes * x7[i] + educationPrimary_or_less * x8[i] +
                         educationUniversity * x9[i] +
                         smokeYes * x10[i] + no2gps_24h * x11[i] + maxwindspeed_24h * x12[i] +
                         precip_24h * x13[i] + sec_noise55_day * x14[i] +
                         access_greenbluespaces_300mbuffYes * x15[i] +
                         age_yrs * x16[i] + tmean_24h * x17[i], sigma^-2)
  }

  # Priors
  alpha ~ dnorm(0, 100^-2)
  genderOtra ~ dnorm(0, 100^-2)
  genderWoman ~ dnorm(0, 100^-2)
  on_a_dietYes ~ dnorm(0, 100^-2)
  alcoholYes ~ dnorm(0, 100^-2)
  drugsYes ~ dnorm(0, 100^-2)
  sickYes ~ dnorm(0, 100^-2)
  other_factorsYes ~ dnorm(0, 100^-2)
  educationPrimary_or_less ~ dnorm(0, 100^-2)
  educationUniversity ~ dnorm(0, 100^-2)
  smokeYes ~ dnorm(0, 100^-2)
  no2gps_24h ~ dnorm(0, 100^-2)
  maxwindspeed_24h ~ dnorm(0, 100^-2)
  precip_24h ~ dnorm(0, 100^-2)
  sec_noise55_day ~ dnorm(0, 100^-2)
  access_greenbluespaces_300mbuffYes ~ dnorm(0, 100^-2)
  age_yrs ~ dnorm(0, 100^-2)
  tmean_24h ~ dnorm(0, 100^-2)
  sigma ~ dunif(0, 10)
}

# Choose the parameters to watch
model_parameters_2a =  c("alpha", "genderOtra","genderWoman","on_a_dietYes","alcoholYes",
                      "drugsYes","sickYes","other_factorsYes", "educationPrimary_or_less",
                      "educationUniversity", "smokeYes","no2gps_24h", "maxwindspeed_24h",
                      "precip_24h","sec_noise55_day","access_greenbluespaces_300mbuffYes",
                      "age_yrs","tmean_24h")

# Set up the data
model_data_2a = list(n = n2_a, y = y2_a,
                  x1 = x2_a[,1], x2 = x2_a[,2], x3 = x2_a[,3],
                  x4 = x2_a[,4], x5 = x2_a[,5], x6 = x2_a[,6],
                  x7 = x2_a[,7], x8 = x2_a[,8], x9 = x2_a[,9],
                  x10 = x2_a[,10], x11 = x2_a[,11], x12 = x2_a[,12],
                  x13 = x2_a[,13], x14 = x2_a[,14], x15 = x2_a[,15],
                  x16 = x2_a[,16], x17 = x2_a[,17])

jags.inits_2a <- function(){
    list('alpha' = stroop_test_performance_lm$coefficients[1],
          'genderOtra' = stroop_test_performance_lm$coefficients[2],
          'genderWoman' = stroop_test_performance_lm$coefficients[3],
          'on_a_dietYes' = stroop_test_performance_lm$coefficients[4],
          'alcoholYes' = stroop_test_performance_lm$coefficients[5],
          'drugsYes' = stroop_test_performance_lm$coefficients[6],
          'sickYes' = stroop_test_performance_lm$coefficients[7],
          'other_factorsYes' = stroop_test_performance_lm$coefficients[8],
          'educationPrimary_or_less' = stroop_test_performance_lm$coefficients[9],
          'educationUniversity' = stroop_test_performance_lm$coefficients[10],
          'smokeYes' = stroop_test_performance_lm$coefficients[11],
          'no2gps_24h' = stroop_test_performance_lm$coefficients[12],
          'maxwindspeed_24h'= stroop_test_performance_lm$coefficients[13],
          'precip_24h' = stroop_test_performance_lm$coefficients[14],
          'sec_noise55_day' = stroop_test_performance_lm$coefficients[15],
          'access_greenbluespaces_300mbuffYes' = stroop_test_performance_lm$coefficients[16],
          'age_yrs' = stroop_test_performance_lm$coefficients[17],
          'tmean_24h' = stroop_test_performance_lm$coefficients[18])
  }

# Run the model
model_run_2a = jags(data = model_data_2a,
                 parameters.to.save = model_parameters_2a,
                 model.file = model_code_2a,
                 inits = jags.inits_2a,
                 n.chains = 1,
                 n.iter = 1000,
                 n.burnin = 200,
                 n.thin = 1)
round(model_run_2a$BUGSoutput$summary,4)

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：JAGS R语言 数据模型

怎么获得积分