{"id":1885,"date":"2025-10-30T07:04:19","date_gmt":"2025-10-29T22:04:19","guid":{"rendered":"https:\/\/info.zanet.biz\/?p=1885"},"modified":"2025-11-01T15:59:01","modified_gmt":"2025-11-01T06:59:01","slug":"stan%e3%82%92%e4%bd%bf%e3%81%86%e3%83%99%e3%82%a4%e3%82%b8%e3%82%a2%e3%83%b3%e3%83%a1%e3%82%bf%e3%82%a2%e3%83%8a%e3%83%aa%e3%82%b7%e3%82%b9%e7%94%a8%e3%81%ae%e3%82%b3%e3%83%bc%e3%83%89%ef%bc%9a","status":"publish","type":"post","link":"https:\/\/info.zanet.biz\/?p=1885","title":{"rendered":"Stan\u3092\u4f7f\u3046\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u306e\u30b3\u30fc\u30c9\uff1a\u4e8c\u5024\u5909\u6570\u30a2\u30a6\u30c8\u30ab\u30e0RR, OR, RD\u7528"},"content":{"rendered":"\n

\u4ee5\u524d\u306e\u6295\u7a3f\u300cStan\u3092\u4f7f\u3046\u30d9\u30a4\u30b8\u30a2\u30f3\u30a2\u30d7\u30ed\u30fc\u30c1\u306e\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\uff1aCopilot\u306b\u805e\u304d\u306a\u304c\u3089<\/a>\u300d\u3067R\u304b\u3089rstan\u3092\u4ecb\u3057\u3066Stan\u3092\u52d5\u304b\u3057\u3001\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3092\u5b9f\u884c\u3059\u308b\u307e\u3067\u306e\u904e\u7a0b\u3092\u7d39\u4ecb\u3057\u307e\u3057\u305f\u3002<\/p>\n\n\n\n

\u4eca\u56de\u3001\u30ea\u30b9\u30af\u6bd4\u3001\u30aa\u30c3\u30ba\u6bd4\u3001\u30ea\u30b9\u30af\u5dee\u3092\u52b9\u679c\u6307\u6a19\u3068\u3059\u308b\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3092\u5b9f\u884c\u3057\u3001Forest plot\u3092\u51fa\u529b\u3059\u308b\u307e\u3067\u306e\u3001\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u4f5c\u6210\u3057\u307e\u3057\u305f\u3002<\/p>\n\n\n\n

\u30c7\u30fc\u30bf\u306f\u3001\u56f3\u306b\u793a\u3059\u3088\u3046\u306a\u5f62\u5f0f\u3067Excel\u3067\u7528\u610f\u3057\u307e\u3059\u3002Risk Ratio\u3001RR\u306e\u30bb\u30eb\u3092Odds Ratio\u3001OR\u3042\u308b\u3044\u306fRisk Difference\u3001 RD\u3068\u66f8\u304d\u63db\u3048\u3066\u52b9\u679c\u6307\u6a19\u306e\u30bf\u30a4\u30d7\u3092\u6307\u5b9a\u3057\u307e\u3059\u3002label\u306e\u305d\u306e\u4ed6\u306e\u9805\u76ee\u3082\u305d\u308c\u305e\u308c\u66f8\u304d\u63db\u3048\u3066\u304f\u3060\u3055\u3044\u3002\u3053\u306e\u4f8b\u3067\u306fUnvaccinated\u306e\u30bb\u30eb\u304c\u5bfe\u7167\u306e\u6cbb\u7642\u3001Vaccinated\u306e\u30bb\u30eb\u304c\u4ecb\u5165\u306e\u6cbb\u7642\u3001TB+\u306e\u30bb\u30eb\u304c\u30a2\u30a6\u30c8\u30ab\u30e0\u306b\u306a\u308a\u307e\u3059\u3002\u52b9\u679c\u6307\u6a19\u306e\u30bf\u30a4\u30d7\u306fRR, OR, RD\u306e\u30e9\u30d9\u30eb\u540d\u3067\u5224\u65ad\u3057\u307e\u3059\u3002<\/p>\n\n\n\n

\"\"<\/figure>\n\n\n\n

\u3053\u306e\u4f8b\u3067\u306f\u3001\u30bb\u30ebB3\u304b\u3089G16\u306e\u7bc4\u56f2\u3092\u9078\u629e\u3057\u3066\uff08\u30c7\u30fc\u30bf\u90e8\u5206\u304c\u4f8b\u3048\u3070\u884c7\u307e\u3067\u3057\u304b\u306a\u3044\u5834\u5408\u306f\u3001\u30bb\u30ebB3\u304b\u3089\u30bb\u30ebG9\u306e\u7bc4\u56f2\u3092\u9078\u629e\u3057\u307e\u3059\uff09\u3001\u4ee5\u4e0b\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092R\u3067\u5b9f\u884c\u3059\u308b\u3068Forest plot\u304c\u51fa\u529b\u3055\u308c\u307e\u3059\u3002\u304b\u306a\u308a\u9577\u3044\u306e\u3067\u30a2\u30b3\u30fc\u30c7\u30a3\u30aa\u30f3\u5f62\u5f0f\u3067\u7573\u3093\u3060\u72b6\u614b\u3067\u3059\u304c\u3001\u53f3\u7aef\u3092\u30af\u30ea\u30c3\u30af\u3059\u308b\u3068\u5c55\u958b\u3055\u308c\u3001\u30b3\u30fc\u30c9\u7528\u306e\u30d5\u30a3\u30fc\u30eb\u30c9\u306b\u8a18\u8ff0\u3057\u3066\u3042\u308a\u307e\u3059\u3002<\/p>\n\n\n\n

\u30b3\u30fc\u30c9\u3092\u3059\u3079\u3066\u3001\u30b3\u30d4\u30fc\u3057\u3066\u3001R\u306eEditor\u306b\u8cbc\u308a\u4ed8\u3051\u307e\u3059\u3002Excel\u3067\u30c7\u30fc\u30bf\u7bc4\u56f2\u3092\u9078\u629e\u3057\u3001\u30b3\u30d4\u30fc\uff08Ctrl+C)\u306e\u5f8c\u3001R\u306b\u623b\u308a\u3001exdat = read.delim(“clipboard”,sep=”\\t”,header=TRUE)\u306e\u884c\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u5b9f\u884c\u3057\u3066\u3001\u30af\u30ea\u30c3\u30d7\u30dc\u30fc\u30c9\u7d4c\u7531\u3067\u30c7\u30fc\u30bf\u3092\u8aad\u307f\u8fbc\u307f\u307e\u3059\u3002\u305d\u306e\u5f8c\u3001\u305d\u308c\u4ee5\u964d\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u3059\u3079\u3066\u5b9f\u884c\u3057\u307e\u3059\u3002<\/p>\n\n\n\n

Stan\u3092\u7528\u3044\u308b\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u30b3\u30fc\u30c9<\/summary>
\n
#####rstan, Stan\u306b\u3088\u308bRR, OR, RD\u306e\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u306e\u30b9\u30af\u30ea\u30d7\u30c8#####\n# \u5fc5\u8981\u306a\u30d1\u30c3\u30b1\u30fc\u30b8\u3092\u30ed\u30fc\u30c9\nlibrary(rstan)\nlibrary(bayesplot)\nlibrary(dplyr)\nlibrary(HDInterval)\nlibrary(forestplot)\nlibrary(metafor)\n\n# Get data via clipboard.\nexdat = read.delim(\"clipboard\",sep=\"\\t\",header=TRUE)\nlabem = exdat[,\"label\"]\nem = labem[6]\nlabe_study = labem[1]\nlabe_cont = labem[2]\nlabe_int = labem[3]\nlabe_outc = labem[4]\nlabe_em = labem[5]\n\n# Format the data.\nexdati = na.omit(exdat)\nif(colnames(exdati)[1]==\"author\"){\nstudy = exdati$author\n}\nif(colnames(exdati)[1]==\"study\"){\nstudy = exdati$study\n}\nnc = exdati$nc\nrc = exdati$rc\nnt = exdati$nt\nrt = exdati$rt\nnum_studies = length(study)\nN = num_studies\nmeta_data <- list(\nN = num_studies,\nstudy = study,\nnc = nc,\nrc = rc,\nnt = nt,\nrt = rt\n)\n\n# Stan\u30e2\u30c7\u30eb\u30b3\u30fc\u30c9\n#RR\nif(em == \"RR\"){\nstan_model_code <- \"data {\n  int<lower=1> N;\n  int nc[N];\n  int rc[N];\n  int nt[N];\n  int rt[N];\n}\n\nparameters {\n  \/\/ \u30ea\u30b9\u30af\u6bd4 (RR) \u7528\n  real logRR[N];\n  real mu_logRR; \/\/ \u30ea\u30b9\u30af\u6bd4\u306e\u7d71\u5408\u5024\uff08\u5bfe\u6570\uff09\n  real<lower=0> sigma_logRR; \/\/ \u30ea\u30b9\u30af\u6bd4\u306e\u7814\u7a76\u9593\u306e\u3070\u3089\u3064\u304d\n\n  real<lower=0, upper=1> p_c[N]; \/\/ \u5bfe\u7167\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u767a\u751f\u78ba\u7387\n}\ntransformed parameters { \/\/ \u5404\u7814\u7a76\u306elogRR, logOR, RD\u3092\u3053\u3053\u306b\u5b9a\u7fa9\n  vector[N] logRR_study; \/\/ \u5404\u7814\u7a76\u306elogRR\n\n  for (i in 1:N) {\n    \/\/ \u3053\u308c\u3089\u306e\u5024\u306fmu_logRR\u306a\u3069\u304b\u3089\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u308b\u3082\u306e\u3067\u3059\u304c\u3001\n    \/\/ \u3053\u3053\u3067\u76f4\u63a5\u89b3\u6e2c\u30c7\u30fc\u30bf\u304b\u3089\u8a08\u7b97\u3055\u308c\u305f\u300c\u70b9\u63a8\u5b9a\u5024\u300d\u3082\u7528\u610f\u3057\u307e\u3059\n    \/\/ \u3053\u308c\u3089\u306fgenerated quantities\u3067\u3057\u304b\u4f7f\u308f\u306a\u3044\u306e\u3067\u3001generated quantities\u306b\u79fb\u3059\u3053\u3068\u3082\u53ef\u80fd\n    \/\/ \u305f\u3060\u3057\u3001model\u30d6\u30ed\u30c3\u30af\u3067\u4f7f\u308f\u308c\u3066\u3044\u308b logRR[i] \u7b49\u3068\u91cd\u8907\u3057\u306a\u3044\u3088\u3046\u306b\u6ce8\u610f\n    \/\/ \u3053\u3053\u3067\u306f\u3001\u8a08\u7b97\u306e\u4fbf\u5b9c\u4e0a\u3001rc,nc,rt,nt \u304b\u3089\u76f4\u63a5\u70b9\u63a8\u5b9a\u5024\u3092\u8a08\u7b97\u3057\u307e\u3059\n    real p_c_obs = (rc[i] + 0.5) \/ (nc[i] + 1.0); \/\/ \u30bc\u30ed\u30a4\u30d9\u30f3\u30c8\u56de\u907f\u306e\u305f\u3081\u306e\u88dc\u6b63\n    real p_t_obs = (rt[i] + 0.5) \/ (nt[i] + 1.0);\n\n    logRR_study[i] = log(p_t_obs) - log(p_c_obs);\n  }\n}\nmodel {\n  \/\/ \u4e8b\u524d\u5206\u5e03\n  \/\/ RR\n  mu_logRR ~ normal(0, 10);\n  sigma_logRR ~ cauchy(0, 2);\n\n  for (i in 1:N) {\n    p_c[i] ~ beta(1, 1); \/\/ \u5bfe\u7167\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u767a\u751f\u78ba\u7387\u306b\u5bfe\u3059\u308b\u4e00\u69d8\u306a\u4e8b\u524d\u5206\u5e03\n\n    \/\/ RR\n    logRR[i] ~ normal(mu_logRR, sigma_logRR);\n    rc[i] ~ binomial(nc[i], p_c[i]);\n    rt[i] ~ binomial(nt[i], fmin(p_c[i] * exp(logRR[i]), 1)); \/\/ RR\u306b\u57fa\u3065\u304d\u4ecb\u5165\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u7387\u3092\u8a08\u7b97\n  }\n}\n\ngenerated quantities {\n  \/\/ \u30ea\u30b9\u30af\u6bd4 (RR)\n  real RR_pooled = exp(mu_logRR);\n  real RR_pred = exp(normal_rng(mu_logRR, sigma_logRR));\n  real tau_squared_logRR = sigma_logRR * sigma_logRR;\n\n  \/\/ \u5404\u7814\u7a76\u306e\u89b3\u6e2c\u3055\u308c\u305f\u52b9\u679c\u91cf\u306e\u6a19\u6e96\u8aa4\u5dee (se_obs) \u3092\u8a08\u7b97\n  vector<lower=0>[N] se_logRR_obs;\n\n  \/\/ \u5404\u52b9\u679c\u6307\u6a19\u306e\u89b3\u6e2c\u8aa4\u5dee\u306e\u5206\u6563 (variance_obs)\n  vector<lower=0>[N] variance_logRR_obs;\n\n  \/\/ \u5178\u578b\u7684\u306a\u7814\u7a76\u5185\u5206\u6563 (S_squared)\n  real S_squared_logRR;\n\n  for (i in 1:N) {\n    \/\/ \u30bc\u30ed\u30a4\u30d9\u30f3\u30c8\/\u7dcf\u6570\u30bc\u30ed\u3092\u907f\u3051\u308b\u305f\u3081\u306b\u30010.5\u306e\u9023\u7d9a\u6027\u88dc\u6b63\u3092\u9069\u7528\n    \/\/ \u983b\u5ea6\u8ad6\u7684\u306a\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3067\u3088\u304f\u4f7f\u308f\u308c\u308b\u624b\u6cd5\u3067\u3059\u3002\n    real current_rc = rc[i] + 0.5;\n    real current_nc = nc[i] + 1.0;\n    real current_rt = rt[i] + 0.5;\n    real current_nt = nt[i] + 1.0;\n\n    \/\/ logRR \u306e\u6a19\u6e96\u8aa4\u5dee\u306e\u8a08\u7b97 (\u30c7\u30eb\u30bf\u6cd5)\n    se_logRR_obs[i] = sqrt(1.0\/current_rc - 1.0\/current_nc + 1.0\/current_rt - 1.0\/current_nt);\n    variance_logRR_obs[i] = se_logRR_obs[i] * se_logRR_obs[i];\n\n  }\n\n  \/\/ \u5178\u578b\u7684\u306a\u7814\u7a76\u5185\u5206\u6563 (S_squared) \u306e\u8a08\u7b97\n  \/\/ \u5404\u7814\u7a76\u306e\u89b3\u6e2c\u8aa4\u5dee\u306e\u5206\u6563\u306e\u5e73\u5747\u5024\n  S_squared_logRR = mean(variance_logRR_obs);\n}\"\n}\n#OR\nif(em == \"OR\"){\nstan_model_code <- \"data {\n  int<lower=1> N;\n  int nc[N];\n  int rc[N];\n  int nt[N];\n  int rt[N];\n}\n\nparameters {\n  \/\/ \u30aa\u30c3\u30ba\u6bd4 (OR) \u7528\n  real logOR[N];\n  real mu_logOR; \/\/ \u30aa\u30c3\u30ba\u6bd4\u306e\u7d71\u5408\u5024\uff08\u5bfe\u6570\uff09\n  real<lower=0> sigma_logOR; \/\/ \u30aa\u30c3\u30ba\u6bd4\u306e\u7814\u7a76\u9593\u306e\u3070\u3089\u3064\u304d\n  real<lower=0, upper=1> p_c[N]; \/\/ \u5bfe\u7167\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u767a\u751f\u78ba\u7387\n}\ntransformed parameters { \/\/ \u5404\u7814\u7a76\u306elogRR, logOR, RD\u3092\u3053\u3053\u306b\u5b9a\u7fa9\n  vector[N] logOR_study; \/\/ \u5404\u7814\u7a76\u306elogOR\n  for (i in 1:N) {\n    \/\/ \u3053\u308c\u3089\u306e\u5024\u306fmu_logRR\u306a\u3069\u304b\u3089\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u308b\u3082\u306e\u3067\u3059\u304c\u3001\n    \/\/ \u3053\u3053\u3067\u76f4\u63a5\u89b3\u6e2c\u30c7\u30fc\u30bf\u304b\u3089\u8a08\u7b97\u3055\u308c\u305f\u300c\u70b9\u63a8\u5b9a\u5024\u300d\u3082\u7528\u610f\u3057\u307e\u3059\n    \/\/ \u3053\u308c\u3089\u306fgenerated quantities\u3067\u3057\u304b\u4f7f\u308f\u306a\u3044\u306e\u3067\u3001generated quantities\u306b\u79fb\u3059\u3053\u3068\u3082\u53ef\u80fd\n    \/\/ \u305f\u3060\u3057\u3001model\u30d6\u30ed\u30c3\u30af\u3067\u4f7f\u308f\u308c\u3066\u3044\u308b logRR[i] \u7b49\u3068\u91cd\u8907\u3057\u306a\u3044\u3088\u3046\u306b\u6ce8\u610f\n    \/\/ \u3053\u3053\u3067\u306f\u3001\u8a08\u7b97\u306e\u4fbf\u5b9c\u4e0a\u3001rc,nc,rt,nt \u304b\u3089\u76f4\u63a5\u70b9\u63a8\u5b9a\u5024\u3092\u8a08\u7b97\u3057\u307e\u3059\n    real p_c_obs = (rc[i] + 0.5) \/ (nc[i] + 1.0); \/\/ \u30bc\u30ed\u30a4\u30d9\u30f3\u30c8\u56de\u907f\u306e\u305f\u3081\u306e\u88dc\u6b63\n    real p_t_obs = (rt[i] + 0.5) \/ (nt[i] + 1.0);\n    logOR_study[i] = log((p_t_obs \/ (1 - p_t_obs))) - log((p_c_obs \/ (1 - p_c_obs)));\n  }\n}\nmodel {\n  \/\/ \u4e8b\u524d\u5206\u5e03\n  \/\/ OR\n  mu_logOR ~ normal(0, 10);\n  sigma_logOR ~ cauchy(0, 2);\n  for (i in 1:N) {\n    p_c[i] ~ beta(1, 1); \/\/ \u5bfe\u7167\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u767a\u751f\u78ba\u7387\u306b\u5bfe\u3059\u308b\u4e00\u69d8\u306a\u4e8b\u524d\u5206\u5e03\n    \/\/ OR\n    logOR[i] ~ normal(mu_logOR, sigma_logOR);\n\trc[i] ~ binomial(nc[i], p_c[i]);\n    \/\/ \u30aa\u30c3\u30ba\u6bd4\u306e\u5b9a\u7fa9 p_t \/ (1 - p_t) = OR * (p_c \/ (1 - p_c)) \u3088\u308a p_t \u3092\u5c0e\u51fa\n    \/\/ logit(p_t) = logOR + logit(p_c)\n    \/\/ p_t = inv_logit(logOR + logit(p_c))\n    rt[i] ~ binomial(nt[i], inv_logit(logOR[i] + logit(p_c[i])));\n  }\n}\ngenerated quantities {\n  \/\/ \u30aa\u30c3\u30ba\u6bd4 (OR)\n  real OR_pooled = exp(mu_logOR);\n  real OR_pred = exp(normal_rng(mu_logOR, sigma_logOR));\n  real tau_squared_logOR = sigma_logOR * sigma_logOR;\n\n  \/\/ \u5404\u7814\u7a76\u306e\u89b3\u6e2c\u3055\u308c\u305f\u52b9\u679c\u91cf\u306e\u6a19\u6e96\u8aa4\u5dee (se_obs) \u3092\u8a08\u7b97\n  vector<lower=0>[N] se_logOR_obs;\n\n  \/\/ \u5404\u52b9\u679c\u6307\u6a19\u306e\u89b3\u6e2c\u8aa4\u5dee\u306e\u5206\u6563 (variance_obs)\n  vector<lower=0>[N] variance_logOR_obs;\n\n  \/\/ \u5178\u578b\u7684\u306a\u7814\u7a76\u5185\u5206\u6563 (S_squared)\n  real S_squared_logOR;\n\n  for (i in 1:N) {\n    \/\/ \u30bc\u30ed\u30a4\u30d9\u30f3\u30c8\/\u7dcf\u6570\u30bc\u30ed\u3092\u907f\u3051\u308b\u305f\u3081\u306b\u30010.5\u306e\u9023\u7d9a\u6027\u88dc\u6b63\u3092\u9069\u7528\n    \/\/ \u983b\u5ea6\u8ad6\u7684\u306a\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3067\u3088\u304f\u4f7f\u308f\u308c\u308b\u624b\u6cd5\u3067\u3059\u3002\n    real current_rc = rc[i] + 0.5;\n    real current_nc = nc[i] + 1.0;\n    real current_rt = rt[i] + 0.5;\n    real current_nt = nt[i] + 1.0;\n\n    \/\/ logOR \u306e\u6a19\u6e96\u8aa4\u5dee\u306e\u8a08\u7b97 (\u30c7\u30eb\u30bf\u6cd5)\n    se_logOR_obs[i] = sqrt(1.0\/current_rc + 1.0\/(current_nc - current_rc) + 1.0\/current_rt + 1.0\/(current_nt - current_rt));\n    variance_logOR_obs[i] = se_logOR_obs[i] * se_logOR_obs[i];\n  }\n  \/\/ \u5178\u578b\u7684\u306a\u7814\u7a76\u5185\u5206\u6563 (S_squared) \u306e\u8a08\u7b97\n  \/\/ \u5404\u7814\u7a76\u306e\u89b3\u6e2c\u8aa4\u5dee\u306e\u5206\u6563\u306e\u5e73\u5747\u5024\n  S_squared_logOR = mean(variance_logOR_obs);\n}\"\n}\n#RD\nif(em == \"RD\"){\nstan_model_code <- \"data {\n  int<lower=1> N;\n  int nc[N];\n  int rc[N];\n  int nt[N];\n  int rt[N];\n}\nparameters {\n  \/\/ \u30ea\u30b9\u30af\u5dee (RD) \u7528\n  real RD[N];\n  real mu_RD; \/\/ \u30ea\u30b9\u30af\u5dee\u306e\u7d71\u5408\u5024\n  real<lower=0> sigma_RD; \/\/ \u30ea\u30b9\u30af\u5dee\u306e\u7814\u7a76\u9593\u306e\u3070\u3089\u3064\u304d\n\n  real<lower=0, upper=1> p_c[N]; \/\/ \u5bfe\u7167\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u767a\u751f\u78ba\u7387\n}\ntransformed parameters { \/\/ \u5404\u7814\u7a76\u306elogRR, logOR, RD\u3092\u3053\u3053\u306b\u5b9a\u7fa9\n  vector[N] RD_study;    \/\/ \u5404\u7814\u7a76\u306eRD\n  for (i in 1:N) {\n    \/\/ \u3053\u308c\u3089\u306e\u5024\u306fmu_logRR\u306a\u3069\u304b\u3089\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u308b\u3082\u306e\u3067\u3059\u304c\u3001\n    \/\/ \u3053\u3053\u3067\u76f4\u63a5\u89b3\u6e2c\u30c7\u30fc\u30bf\u304b\u3089\u8a08\u7b97\u3055\u308c\u305f\u300c\u70b9\u63a8\u5b9a\u5024\u300d\u3082\u7528\u610f\u3057\u307e\u3059\n    \/\/ \u3053\u308c\u3089\u306fgenerated quantities\u3067\u3057\u304b\u4f7f\u308f\u306a\u3044\u306e\u3067\u3001generated quantities\u306b\u79fb\u3059\u3053\u3068\u3082\u53ef\u80fd\n    \/\/ \u305f\u3060\u3057\u3001model\u30d6\u30ed\u30c3\u30af\u3067\u4f7f\u308f\u308c\u3066\u3044\u308b logRR[i] \u7b49\u3068\u91cd\u8907\u3057\u306a\u3044\u3088\u3046\u306b\u6ce8\u610f\n    \/\/ \u3053\u3053\u3067\u306f\u3001\u8a08\u7b97\u306e\u4fbf\u5b9c\u4e0a\u3001rc,nc,rt,nt \u304b\u3089\u76f4\u63a5\u70b9\u63a8\u5b9a\u5024\u3092\u8a08\u7b97\u3057\u307e\u3059\n    real p_c_obs = (rc[i] + 0.5) \/ (nc[i] + 1.0); \/\/ \u30bc\u30ed\u30a4\u30d9\u30f3\u30c8\u56de\u907f\u306e\u305f\u3081\u306e\u88dc\u6b63\n    real p_t_obs = (rt[i] + 0.5) \/ (nt[i] + 1.0);\n    RD_study[i] = p_t_obs - p_c_obs;\n  }\n}\nmodel {\n  \/\/ \u4e8b\u524d\u5206\u5e03\n  \/\/ RD\n  mu_RD ~ normal(0, 10);\n  sigma_RD ~ cauchy(0, 2);\n  for (i in 1:N) {\n    p_c[i] ~ beta(1, 1); \/\/ \u5bfe\u7167\u7fa4\u306e\u30a4\u30d9\u30f3\u30c8\u767a\u751f\u78ba\u7387\u306b\u5bfe\u3059\u308b\u4e00\u69d8\u306a\u4e8b\u524d\u5206\u5e03\n\n    \/\/ RD\n    RD[i] ~ normal(mu_RD, sigma_RD);\n\trc[i] ~ binomial(nc[i], p_c[i]);\n    \/\/ \u30ea\u30b9\u30af\u5dee\u306e\u5b9a\u7fa9 p_t - p_c = RD \u3088\u308a p_t \u3092\u5c0e\u51fa\n    \/\/ p_t = p_c + RD\n    rt[i] ~ binomial(nt[i], fmin(fmax(p_c[i] + RD[i], 0), 1)); \/\/ 0\u304b\u30891\u306e\u7bc4\u56f2\u306b\u5236\u9650\n  }\n}\ngenerated quantities {\n  \/\/ \u30ea\u30b9\u30af\u5dee (RD)\n  real RD_pooled = mu_RD;\n  real RD_pred = normal_rng(mu_RD, sigma_RD);\n  real tau_squared_RD = sigma_RD * sigma_RD;\n  \n  \/\/ \u5404\u7814\u7a76\u306e\u89b3\u6e2c\u3055\u308c\u305f\u52b9\u679c\u91cf\u306e\u6a19\u6e96\u8aa4\u5dee (se_obs) \u3092\u8a08\u7b97\n  vector<lower=0>[N] se_RD_obs;\n\n  \/\/ \u5404\u52b9\u679c\u6307\u6a19\u306e\u89b3\u6e2c\u8aa4\u5dee\u306e\u5206\u6563 (variance_obs)\n  vector<lower=0>[N] variance_RD_obs;\n\n  \/\/ \u5178\u578b\u7684\u306a\u7814\u7a76\u5185\u5206\u6563 (S_squared)\n  real S_squared_RD;\n  for (i in 1:N) {\n    \/\/ \u30bc\u30ed\u30a4\u30d9\u30f3\u30c8\/\u7dcf\u6570\u30bc\u30ed\u3092\u907f\u3051\u308b\u305f\u3081\u306b\u30010.5\u306e\u9023\u7d9a\u6027\u88dc\u6b63\u3092\u9069\u7528\n    \/\/ \u983b\u5ea6\u8ad6\u7684\u306a\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3067\u3088\u304f\u4f7f\u308f\u308c\u308b\u624b\u6cd5\u3067\u3059\u3002\n    real current_rc = rc[i] + 0.5;\n    real current_nc = nc[i] + 1.0;\n    real current_rt = rt[i] + 0.5;\n    real current_nt = nt[i] + 1.0;\n\n    \/\/ RD \u306e\u6a19\u6e96\u8aa4\u5dee\u306e\u8a08\u7b97 (\u30c7\u30eb\u30bf\u6cd5)\n    \/\/ p_c = rc\/nc, p_t = rt\/nt\n    \/\/ var(RD) = var(p_t - p_c) = var(p_t) + var(p_c) (\u72ec\u7acb\u3068\u4eee\u5b9a)\n    \/\/ var(p) = p*(1-p)\/n\n    real p_c_calc = current_rc \/ current_nc;\n    real p_t_calc = current_rt \/ current_nt;\n    se_RD_obs[i] = sqrt(p_c_calc * (1.0 - p_c_calc) \/ current_nc + p_t_calc * (1.0 - p_t_calc) \/ current_nt);\n    variance_RD_obs[i] = se_RD_obs[i] * se_RD_obs[i];\n  }\n  \/\/ \u5178\u578b\u7684\u306a\u7814\u7a76\u5185\u5206\u6563 (S_squared) \u306e\u8a08\u7b97\n  \/\/ \u5404\u7814\u7a76\u306e\u89b3\u6e2c\u8aa4\u5dee\u306e\u5206\u6563\u306e\u5e73\u5747\u5024\n  S_squared_RD = mean(variance_RD_obs);\n}\"\n}\n# \u521d\u671f\u5024\u306e\u8a2d\u5b9a\nN <- meta_data$N\ninit_values <- function() {\n  list(\n    mu_logRR = 0, sigma_logRR = 1, logRR = rep(0, N),\n    mu_logOR = 0, sigma_logOR = 1, logOR = rep(0, N),\n    mu_RD = 0, sigma_RD = 1, RD = rep(0, N),\n    p_c = rep(0.1, N)\n  )\n}\n# Stan\u30e2\u30c7\u30eb\u306e\u30b3\u30f3\u30d1\u30a4\u30eb\u3068\u5b9f\u884c\n# \u30a6\u30a9\u30fc\u30e0\u30a2\u30c3\u30d710000\u56de\u3001\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u306f20000\u56de\uff58\uff14\u30c1\u30a7\u30a4\u30f3\u306780000\u500b\u306b\u8a2d\u5b9a\u3002\u9069\u5b9c\u5909\u66f4\u3059\u308b\u3002\nfit <- stan(model_code = stan_model_code, data = meta_data, init = init_values, iter = 30000, warmup=10000, chains = 4, seed = 123,control = list(adapt_delta = 0.98))\n#fit <- stan(model_code = stan_model_code, data = meta_data, init = init_values, iter = 30000, warmup=10000, chains = 4, seed = 123)\n\n# \u63a8\u5b9a\u7d50\u679c\u306e\u62bd\u51fa\nposterior_samples <- extract(fit)\n\n# --- \u7d50\u679c\u306e\u8868\u793a ---\n## \u30ea\u30b9\u30af\u6bd4 (RR)\nif(em == \"RR\"){\ncat(\"--- \u30ea\u30b9\u30af\u6bd4 (RR) \u306e\u7d50\u679c ---\\n\")\nmu_logRR_samples <- posterior_samples$mu_logRR\nrisk_ratio <- exp(mean(mu_logRR_samples))\nconf_int_rr <- exp(quantile(mu_logRR_samples, probs = c(0.025, 0.975)))\npi_rr = quantile(posterior_samples$RR_pred, probs = c(0.025, 0.975))\t#PI of RR\np_val_rr = round(2*min(mean(mu_logRR_samples >0), mean(mu_logRR_samples <0)), 6)\t#Bayesian p-value for summary logRR\nif(mean(mu_logRR_samples)>0){\nprob_direct_rr = round(mean(mu_logRR_samples > 0),6)\nlabel = \"p(RR>1)=\"\n}else{\nprob_direct_rr = round(mean(mu_logRR_samples < 0),6)\nlabel = \"p(RR<1)=\"\n}\ncat(\"\u63a8\u5b9a\u3055\u308c\u305f\u30ea\u30b9\u30af\u6bd4:\", risk_ratio, \"\\n\")\ncat(\"95% \u78ba\u4fe1\u533a\u9593:\", conf_int_rr[1], \"-\", conf_int_rr[2],\" P=\",p_val_rr,label,prob_direct_rr,\"\\n\\n\")\ncat(\"95% \u63a8\u6e2c\u533a\u9593\uff1a\", pi_rr[1], \"-\", pi_rr[2])\n}\nif(em == \"OR\"){\n## \u30aa\u30c3\u30ba\u6bd4 (OR)\ncat(\"--- \u30aa\u30c3\u30ba\u6bd4 (OR) \u306e\u7d50\u679c ---\\n\")\nmu_logOR_samples <- posterior_samples$mu_logOR\nodds_ratio <- exp(mean(mu_logOR_samples))\nconf_int_or <- exp(quantile(mu_logOR_samples, probs = c(0.025, 0.975)))\npi_or = quantile(posterior_samples$OR_pred, probs = c(0.025, 0.975))\t#PI of OR\np_val_or = round(2*min(mean(mu_logOR_samples >0), mean(mu_logOR_samples <0)), 6)\t#Bayesian p-value for summary logOR\nif(mean(mu_logOR_samples)>0){\nprob_direct_or = round(mean(mu_logOR_samples > 0),6)\nlabel = \"p(OR>1)=\"\n}else{\nprob_direct_or = round(mean(mu_logOR_samples < 0),6)\nlabel = \"p(OR<1)=\"\n}\ncat(\"\u63a8\u5b9a\u3055\u308c\u305f\u30aa\u30c3\u30ba\u6bd4:\", odds_ratio, \"\\n\")\ncat(\"95% \u78ba\u4fe1\u533a\u9593:\", conf_int_or[1], \"-\", conf_int_or[2],\" P=\",p_val_or,label,prob_direct_or, \"\\n\\n\")\ncat(\"95% \u63a8\u6e2c\u533a\u9593\uff1a\", pi_or[1], \"-\", pi_or[2])\n}\nif(em == \"RD\"){\n## \u30ea\u30b9\u30af\u5dee (RD)\ncat(\"--- \u30ea\u30b9\u30af\u5dee (RD) \u306e\u7d50\u679c ---\\n\")\nmu_RD_samples <- posterior_samples$mu_RD\nrisk_difference <- mean(mu_RD_samples)\nconf_int_rd <- quantile(mu_RD_samples, probs = c(0.025, 0.975))\npi_rd = quantile(posterior_samples$RD_pred, probs = c(0.025, 0.975))\t#PI of RD\np_val_rd = round(2*min(mean(mu_RD_samples >0), mean(mu_RD_samples <0)), 6)\t#Bayesian p-value for summary RD\nif(mean(mu_RD_samples)>0){\nprob_direct_rd = round(mean(mu_RD_samples > 0),6)\nlabel = \"p(RD>0)=\"\n}else{\nprob_direct_rd = round(mean(mu_RD_samples < 0),6)\nlabel = \"p(RD<0)=\"\n}\ncat(\"\u63a8\u5b9a\u3055\u308c\u305f\u30ea\u30b9\u30af\u5dee:\", risk_difference, \"\\n\")\ncat(\"95% \u78ba\u4fe1\u533a\u9593:\", conf_int_rd[1], \"-\", conf_int_rd[2],\" P=\",p_val_rd,label,prob_direct_rd, \"\\n\")\ncat(\"95% \u63a8\u6e2c\u533a\u9593\uff1a\", pi_rd[1], \"-\", pi_rd[2])\n}\n#######\n#\u5404\u7814\u7a76\u306e\u52b9\u679c\u6307\u6a19\u306895%\u78ba\u4fe1\u533a\u9593\u306e\u5024\nif(em == \"RR\"){\n# \u30ea\u30b9\u30af\u6bd4\nrr=rep(0,N)\nrr_lw=rep(0,N)\nrr_up=rep(0,N)\nfor(i in 1:N){\nrr[i] = exp(mean(posterior_samples$logRR[,i]))\nrr_lw[i] = exp(quantile(posterior_samples$logRR[,i],probs=0.025))\nrr_up[i] = exp(quantile(posterior_samples$logRR[,i],probs=0.975))\n}\n}\nif(em == \"OR\"){\n#\u30aa\u30c3\u30ba\u6bd4\nor=rep(0,N)\nor_lw=rep(0,N)\nor_up=rep(0,N)\nfor(i in 1:N){\nor[i] = exp(mean(posterior_samples$logOR[,i]))\nor_lw[i] = exp(quantile(posterior_samples$logOR[,i],probs=0.025))\nor_up[i] = exp(quantile(posterior_samples$logOR[,i],probs=0.975))\n}\n}\nif(em == \"RD\"){\n#\u30ea\u30b9\u30af\u5dee\nrd=rep(0,N)\nrd_lw=rep(0,N)\nrd_up=rep(0,N)\nfor(i in 1:N){\nrd[i] = mean(posterior_samples$RD[,i])\nrd_lw[i] = quantile(posterior_samples$RD[,i],probs=0.025)\nrd_up[i] = quantile(posterior_samples$RD[,i],probs=0.975)\n}\n}\n#####\n#####I-squared:\nif(em == \"RR\"){\n# --- \u30ea\u30b9\u30af\u6bd4 (logRR) \u306e I-squared \u8a08\u7b97 ---\ntau_squared_logRR_samples <- posterior_samples$tau_squared_logRR\nS_squared_logRR_from_stan <- posterior_samples$S_squared_logRR[1] # S_squared\u306f\u5404\u30a4\u30c6\u30ec\u30fc\u30b7\u30e7\u30f3\u3067\u540c\u3058\u5024\u306a\u306e\u3067\u3001\u6700\u521d\u306e\u8981\u7d20\u3067OK\n\n# \u5404\u30a4\u30c6\u30ec\u30fc\u30b7\u30e7\u30f3\u3067I-squared\u3092\u8a08\u7b97\nI_squared_logRR_samples <- (tau_squared_logRR_samples \/ (tau_squared_logRR_samples + S_squared_logRR_from_stan)) * 100\n\n# I-squared\u306e\u70b9\u63a8\u5b9a\u5024\uff08\u5e73\u5747\u5024\uff09\u306895%\u78ba\u4fe1\u533a\u9593\u3092\u8a08\u7b97\n#mean_I_squared_logRR <- mean(I_squared_logRR_samples)\n#ci_I_squared_logRR <- quantile(I_squared_logRR_samples, probs = c(0.025, 0.975))\n#mode and HDI (high density interval) RR\ndens=density(I_squared_logRR_samples)\nmean_I_squared_logRR=dens$x[which.max(dens$y)]\t#mode\nci_I_squared_logRR=hdi(I_squared_logRR_samples)\t#High Density Interval (HDI) lower and upper.\n\ncat(paste0(\"--- I-squared for Log Relative Risk ---\\n\"))\ncat(paste0(\"Estimated I-squared (Mode): \", round(mean_I_squared_logRR, 2), \"%\\n\"))\ncat(paste0(\"95% High Density Interval for I-squared: \", round(ci_I_squared_logRR[1], 2), \"% - \", round(ci_I_squared_logRR[2], 2), \"%\\n\"))\n}\nif(em == \"OR\"){\n# --- \u30aa\u30c3\u30ba\u6bd4 (logOR)\u3064\u3044\u3066\u3082\u540c\u69d8\u306b\u8a08\u7b97 ---\ntau_squared_logOR_samples <- posterior_samples$tau_squared_logOR\nS_squared_logOR_from_stan <- posterior_samples$S_squared_logOR[1] # S_squared\u306f\u5404\u30a4\u30c6\u30ec\u30fc\u30b7\u30e7\u30f3\u3067\u540c\u3058\u5024\u306a\u306e\u3067\u3001\u6700\u521d\u306e\u8981\u7d20\u3067OK\n# ... \u540c\u69d8\u306e\u8a08\u7b97 ...\n# \u5404\u30a4\u30c6\u30ec\u30fc\u30b7\u30e7\u30f3\u3067I-squared\u3092\u8a08\u7b97\nI_squared_logOR_samples <- (tau_squared_logOR_samples \/ (tau_squared_logOR_samples + S_squared_logOR_from_stan)) * 100\n\n# I-squared\u306e\u70b9\u63a8\u5b9a\u5024\uff08\u5e73\u5747\u5024\uff09\u306895%\u78ba\u4fe1\u533a\u9593\u3092\u8a08\u7b97\n#mean_I_squared_logOR <- mean(I_squared_logOR_samples)\n#ci_I_squared_logOR <- quantile(I_squared_logOR_samples, probs = c(0.025, 0.975))\n#mode and HDI (High Density Interval) OR\ndens=density(I_squared_logOR_samples)\nmean_I_squared_logOR=dens$x[which.max(dens$y)]\t#mode\nci_I_squared_logOR=hdi(I_squared_logOR_samples)\t#High Density Interval (HDI) lower and upper.\n\ncat(paste0(\"--- I-squared for Log Odds Ratio ---\\n\"))\ncat(paste0(\"Estimated I-squared (Mode): \", round(mean_I_squared_logOR, 2), \"%\\n\"))\ncat(paste0(\"95% High Density Interval for I-squared: \", round(ci_I_squared_logOR[1], 2), \"% - \", round(ci_I_squared_logOR[2], 2), \"%\\n\"))\n}\nif(em == \"RD\"){\n# --- \u30ea\u30b9\u30af\u5dee (RD) \u306b\u3064\u3044\u3066\u3082\u540c\u69d8\u306b\u8a08\u7b97 ---\ntau_squared_RD_samples <- posterior_samples$tau_squared_RD\nS_squared_RD_from_stan <- posterior_samples$S_squared_RD[1] # S_squared\u306f\u5404\u30a4\u30c6\u30ec\u30fc\u30b7\u30e7\u30f3\u3067\u540c\u3058\u5024\u306a\u306e\u3067\u3001\u6700\u521d\u306e\u8981\u7d20\u3067OK\n# ... \u540c\u69d8\u306e\u8a08\u7b97 ...\n# \u5404\u30a4\u30c6\u30ec\u30fc\u30b7\u30e7\u30f3\u3067I-squared\u3092\u8a08\u7b97\nI_squared_RD_samples <- (tau_squared_RD_samples \/ (tau_squared_RD_samples + S_squared_RD_from_stan)) * 100\n\n# I-squared\u306e\u70b9\u63a8\u5b9a\u5024\uff08\u5e73\u5747\u5024\uff09\u306895%\u78ba\u4fe1\u533a\u9593\u3092\u8a08\u7b97\n#mean_I_squared_RD <- mean(I_squared_RD_samples)\n#ci_I_squared_RD <- quantile(I_squared_RD_samples, probs = c(0.025, 0.975))\n#mode and HDI (High Density Interval) RD\ndens=density(I_squared_RD_samples)\nmean_I_squared_RD=dens$x[which.max(dens$y)]\t#mode\nci_I_squared_RD=hdi(I_squared_RD_samples)\t#High Density Interval (HDI) lower and upper.\n\ncat(paste0(\"--- I-squared for Risk Difference ---\\n\"))\ncat(paste0(\"Estimated I-squared (Mode): \", round(mean_I_squared_RD, 2), \"%\\n\"))\ncat(paste0(\"95% High Density Interval for I-squared: \", round(ci_I_squared_RD[1], 2), \"% - \", round(ci_I_squared_RD[2], 2), \"%\\n\"))\n}\n#####\n####MCMC diagnostic####\n#MCMC trace plot: logRR\n#dev.new();mcmc_dens(as.matrix(fit, pars = c(\"logRR\",\"tau_squared_logRR\")))\n#dev.new();mcmc_trace(fit, pars = c(\"mu_logRR\", \"sigma_logRR\"))\n#Pairs plot\n#dev.new();mcmc_pairs(fit, pars = c(\"mu_logRR\", \"sigma_logRR\"))\n\n#### Forest plot #####\n#Risk Ratio (RR)\nif(em == \"RR\"){\n# \u91cd\u307f\u3092\u683c\u7d0d\u3059\u308b\u30d9\u30af\u30c8\u30eb\u3092\u521d\u671f\u5316\nweights <- numeric(N)\nweight_percentages <- numeric(N)\n\nfor (i in 1:N) {\n  # i\u756a\u76ee\u306e\u7814\u7a76\u306elogRR\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u5024\u3092\u53d6\u5f97\n  log_RR_i_samples <- posterior_samples$logRR[, i]\n  # \u305d\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u305f\u5024\u306e\u5206\u6563\u3092\u8a08\u7b97\n  # \u3053\u308c\u3092\u300c\u5b9f\u52b9\u7684\u306a\u9006\u5206\u6563\u300d\u3068\u8003\u3048\u308b\n  variance_of_logRR_i <- var(log_RR_i_samples)\n  # \u91cd\u307f\u306f\u5206\u6563\u306e\u9006\u6570\n  weights[i] <- 1 \/ variance_of_logRR_i\n}\n# \u5168\u91cd\u307f\u306e\u5408\u8a08\ntotal_weight <- sum(weights)\n# \u5404\u7814\u7a76\u306e\u91cd\u307f\u306e\u30d1\u30fc\u30bb\u30f3\u30c6\u30fc\u30b8\nweight_percentages <- (weights \/ total_weight) * 100\n# \u7d50\u679c\u306e\u8868\u793a\nresults_weights_posterior_var <- data.frame(\n  Study = 1:N,\n  Effective_Weight = weights,\n  Effective_Weight_Percentage = weight_percentages\n)\nprint(results_weights_posterior_var)\nk=N\nwpc = format(round(weight_percentages,digits=1), nsmall=1)\nwp = weight_percentages\/100\n#Forest plot box sizes on weihts\nwbox=c(NA,NA,(k\/4)*sqrt(wp)\/sum(sqrt(wp)),0.5,0,NA)\t\n#wbox=c(NA,NA,(k\/5)*sqrt(wp)\/sum(sqrt(wp)),0.5,0,NA)\n###### tau-squared\n# logRR \u306e tau^2\ntau_squared_logRR_samples <- posterior_samples$tau_squared_logRR\n#mean_tau_squared_logRR <- mean(tau_squared_logRR_samples)\n#ci_tau_squared_logRR <- quantile(tau_squared_logRR_samples, probs = c(0.025, 0.975))\n#mode and HDI (High Density Interval) logRR tau-squared\ndens=density(tau_squared_logRR_samples)\nmean_tau_squared_logRR=dens$x[which.max(dens$y)]\t#mode\nci_tau_squared_logRR=hdi(tau_squared_logRR_samples)\t#High Density Interval (HDI) lower and upper.\n\nprint(paste(\"Estimated tau-squared (logRR, mode):\", mean_tau_squared_logRR))\nprint(paste(\"95% High Density Interval for tau-squared (logRR):\", ci_tau_squared_logRR[1], \"-\", ci_tau_squared_logRR[2]))\n\n#Forest plot by forestplot\n#setting fs for cex\nfs=1\nif(k>20){fs=round((1-0.02*(k-20)),digits=1)}\n\nm=c(NA,NA,rr,risk_ratio,risk_ratio,NA)\nlw=c(NA,NA,rr_lw,conf_int_rr[1],pi_rr[1],NA)\nup=c(NA,NA,rr_up,conf_int_rr[2],pi_rr[2],NA)\n\nhete1=\"\"\nhete2=paste(\"I2=\",round(mean_I_squared_logRR, 2),\"%\",sep=\"\")\n\nhete3=paste(\"tau2=\",format(round(mean_tau_squared_logRR,digits=4),nsmall=4),sep=\"\")\n\n#hete2=paste(\"I2=\",round(mean_I_squared_logRR, 2), \"%\\n\",\" (\",round(ci_I_squared_logRR[1], 2),\" ~ \",\n#round(ci_I_squared_logRR[2], 2),\")\",sep=\"\")\n\n#hete3=paste(\"tau2=\",format(round(mean_tau_squared_logRR,digits=4),nsmall=4), \"\\n\",\n#\" (\",format(round(ci_tau_squared_logRR[1],digits=4),nsmall=4),\" ~ \",\n#format(round(ci_tau_squared_logRR[2],digits=4),nsmall=4),\")  \",sep=\"\")\n\nhete4=\"\"\nhete5=paste(\"p=\",p_val_rr, sep=\"\")\nhete6=paste(label,prob_direct_rr,sep=\"\")\n\nau=study\nsl=c(NA,toString(labe_study),as.vector(au),\"Summary Estimate\",\"Prediction Interval\",NA)\n\nncl=c(labe_cont,\"Number\",nc,NA,NA,hete1)\nrcl=c(NA,labe_outc,rc,NA,NA,hete2)\n\nntl=c(labe_int,\"Number\",nt,NA,NA,hete3)\nrtl=c(NA,labe_outc,rt,NA,NA,hete4)\n\nspac=c(\"    \",NA,rep(NA,k),NA,NA,NA)\nml=c(NA,labe_em,format(round(rr,digits=3),nsmall=3),format(round(risk_ratio,digits=3),nsmall=3),NA,hete5)\nll=c(NA,\"95% CI lower\",format(round(rr_lw,digits=3),nsmall=3),format(round(conf_int_rr[1],digits=3),nsmall=3),format(round(pi_rr[1],digits=3),nsmall=3),hete6)\nul=c(NA,\"95% CI upper\",format(round(rr_up,digits=3),nsmall=3),format(round(conf_int_rr[2],digits=3),nsmall=3),format(round(pi_rr[2],digits=3),nsmall=3),NA)\nwpcl=c(NA,\"Weight(%)\",wpc,100,NA,NA)\nll=as.vector(ll)\nul=as.vector(ul)\nsum=c(TRUE,TRUE,rep(FALSE,k),TRUE,FALSE,FALSE)\nzerov=1\nlabeltext=cbind(sl,ntl,rtl,ncl,rcl,spac,ml,ll,ul,wpcl)\nhlines=list(\"3\"=gpar(lwd=1,columns=1:11,col=\"grey\"))\ndev.new()\nplot(forestplot(labeltext,mean=m,lower=lw,upper=up,is.summary=sum,graph.pos=7,\nzero=zerov,hrzl_lines=hlines,xlab=toString(labe_em),txt_gp=fpTxtGp(ticks=gpar(cex=fs),\nxlab=gpar(cex=fs),cex=fs),xticks.digits=2,vertices=TRUE,graphwidth=unit(50,\"mm\"),colgap=unit(3,\"mm\"),\nboxsize=wbox,\nlineheight=\"auto\",xlog=TRUE,new_page=FALSE))\n}\n#clip=c(clipl,clipu),\n##############################################\n#Odds Ratio (OR)\nif(em == \"OR\"){\n# \u91cd\u307f\u3092\u683c\u7d0d\u3059\u308b\u30d9\u30af\u30c8\u30eb\u3092\u521d\u671f\u5316\nweights <- numeric(N)\nweight_percentages <- numeric(N)\n\nfor (i in 1:N) {\n  # i\u756a\u76ee\u306e\u7814\u7a76\u306elogRR\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u5024\u3092\u53d6\u5f97\n  log_OR_i_samples <- posterior_samples$logOR[, i]\n  # \u305d\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u305f\u5024\u306e\u5206\u6563\u3092\u8a08\u7b97\n  # \u3053\u308c\u3092\u300c\u5b9f\u52b9\u7684\u306a\u9006\u5206\u6563\u300d\u3068\u8003\u3048\u308b\n  variance_of_logOR_i <- var(log_OR_i_samples)\n  # \u91cd\u307f\u306f\u5206\u6563\u306e\u9006\u6570\n  weights[i] <- 1 \/ variance_of_logOR_i\n}\n# \u5168\u91cd\u307f\u306e\u5408\u8a08\ntotal_weight <- sum(weights)\n# \u5404\u7814\u7a76\u306e\u91cd\u307f\u306e\u30d1\u30fc\u30bb\u30f3\u30c6\u30fc\u30b8\nweight_percentages <- (weights \/ total_weight) * 100\n# \u7d50\u679c\u306e\u8868\u793a\nresults_weights_posterior_var <- data.frame(\n  Study = 1:N,\n  Effective_Weight = weights,\n  Effective_Weight_Percentage = weight_percentages\n)\nprint(results_weights_posterior_var)\nk=N\nwpc = format(round(weight_percentages,digits=1), nsmall=1)\nwp = weight_percentages\/100\n#Forest plot box sizes on weihts\nwbox=c(NA,NA,(k\/4)*sqrt(wp)\/sum(sqrt(wp)),0.5,0,NA)\t\n#wbox=c(NA,NA,(k\/5)*sqrt(wp)\/sum(sqrt(wp)),0.5,0,NA)\n###### tau-squared\n# logOR \u306e tau^2\ntau_squared_logOR_samples <- posterior_samples$tau_squared_logOR\n#mean_tau_squared_logOR <- mean(tau_squared_logOR_samples)\n#ci_tau_squared_logOR <- quantile(tau_squared_logOR_samples, probs = c(0.025, 0.975))\n#mode and HDI (High Density Interval) logOR tau-squared\ndens=density(tau_squared_logOR_samples)\nmean_tau_squared_logOR=dens$x[which.max(dens$y)]\t#mode\nci_tau_squared_logOR=hdi(tau_squared_logOR_samples)\t#High Density Interval (HDI) lower and upper.\n\nprint(paste(\"Estimated tau-squared (logOR, mode):\", mean_tau_squared_logOR))\nprint(paste(\"95% High Density Interval for tau-squared (logOR):\", ci_tau_squared_logOR[1], \"-\", ci_tau_squared_logOR[2]))\n\n#Forest plot by forestplot\n#setting fs for cex\nfs=1\nif(k>20){fs=round((1-0.02*(k-20)),digits=1)}\n\nm=c(NA,NA,or,odds_ratio,odds_ratio,NA)\nlw=c(NA,NA,or_lw,conf_int_or[1],pi_or[1],NA)\nup=c(NA,NA,or_up,conf_int_or[2],pi_or[2],NA)\n\nhete1=\"\"\nhete2=paste(\"I2=\",round(mean_I_squared_logOR, 2),\"%\",sep=\"\")\n\nhete3=paste(\"tau2=\",format(round(mean_tau_squared_logOR,digits=4),nsmall=4),sep=\"\")\n\n#hete2=paste(\"I2=\",round(mean_I_squared_logOR, 2), \"%\\n\",\" (\",round(ci_I_squared_logOR[1], 2),\" ~ \",\n#round(ci_I_squared_logOR[2], 2),\")\",sep=\"\")\n\n#hete3=paste(\"tau2=\",format(round(mean_tau_squared_logOR,digits=4),nsmall=4), \"\\n\",\n#\" (\",format(round(ci_tau_squared_logOR[1],digits=4),nsmall=4),\" ~ \",\n#format(round(ci_tau_squared_logOR[2],digits=4),nsmall=4),\")  \",sep=\"\")\n\nhete4=\"\"\nhete5=paste(\"p=\",p_val_or, sep=\"\")\nhete6=paste(label,prob_direct_or,sep=\"\")\n\nau=study\nsl=c(NA,toString(labe_study),as.vector(au),\"Summary Estimate\",\"Prediction Interval\",NA)\n\nncl=c(labe_cont,\"Number\",nc,NA,NA,hete1)\nrcl=c(NA,labe_outc,rc,NA,NA,hete2)\n\nntl=c(labe_int,\"Number\",nt,NA,NA,hete3)\nrtl=c(NA,labe_outc,rt,NA,NA,hete4)\n\nspac=c(\"    \",NA,rep(NA,k),NA,NA,NA)\nml=c(NA,labe_em,format(round(or,digits=3),nsmall=3),format(round(odds_ratio,digits=3),nsmall=3),NA,hete5)\nll=c(NA,\"95% CI lower\",format(round(or_lw,digits=3),nsmall=3),format(round(conf_int_or[1],digits=3),nsmall=3),format(round(pi_or[1],digits=3),nsmall=3),hete6)\nul=c(NA,\"95% CI upper\",format(round(or_up,digits=3),nsmall=3),format(round(conf_int_or[2],digits=3),nsmall=3),format(round(pi_or[2],digits=3),nsmall=3),NA)\nwpcl=c(NA,\"Weight(%)\",wpc,100,NA,NA)\nll=as.vector(ll)\nul=as.vector(ul)\nsum=c(TRUE,TRUE,rep(FALSE,k),TRUE,FALSE,FALSE)\nzerov=1\nlabeltext=cbind(sl,ntl,rtl,ncl,rcl,spac,ml,ll,ul,wpcl)\nhlines=list(\"3\"=gpar(lwd=1,columns=1:11,col=\"grey\"))\ndev.new()\nplot(forestplot(labeltext,mean=m,lower=lw,upper=up,is.summary=sum,graph.pos=7,\nzero=zerov,hrzl_lines=hlines,xlab=toString(labe_em),txt_gp=fpTxtGp(ticks=gpar(cex=fs),\nxlab=gpar(cex=fs),cex=fs),xticks.digits=2,vertices=TRUE,graphwidth=unit(50,\"mm\"),colgap=unit(3,\"mm\"),\nboxsize=wbox,\nlineheight=\"auto\",xlog=TRUE,new_page=FALSE))\n}\n\n#############################################\n#####Risk Difference\nif(em == \"RD\"){\n# \u91cd\u307f\u3092\u683c\u7d0d\u3059\u308b\u30d9\u30af\u30c8\u30eb\u3092\u521d\u671f\u5316\nweights <- numeric(N)\nweight_percentages <- numeric(N)\n\nfor (i in 1:N) {\n  # i\u756a\u76ee\u306e\u7814\u7a76\u306elogRR\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u5024\u3092\u53d6\u5f97\n  RD_i_samples <- posterior_samples$RD[, i]\n  # \u305d\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u305f\u5024\u306e\u5206\u6563\u3092\u8a08\u7b97\n  # \u3053\u308c\u3092\u300c\u5b9f\u52b9\u7684\u306a\u9006\u5206\u6563\u300d\u3068\u8003\u3048\u308b\n  variance_of_RD_i <- var(RD_i_samples)\n  # \u91cd\u307f\u306f\u5206\u6563\u306e\u9006\u6570\n  weights[i] <- 1 \/ variance_of_RD_i\n}\n# \u5168\u91cd\u307f\u306e\u5408\u8a08\ntotal_weight <- sum(weights)\n# \u5404\u7814\u7a76\u306e\u91cd\u307f\u306e\u30d1\u30fc\u30bb\u30f3\u30c6\u30fc\u30b8\nweight_percentages <- (weights \/ total_weight) * 100\n# \u7d50\u679c\u306e\u8868\u793a\nresults_weights_posterior_var <- data.frame(\n  Study = 1:N,\n  Effective_Weight = weights,\n  Effective_Weight_Percentage = weight_percentages\n)\nprint(results_weights_posterior_var)\nk=N\nwpc = format(round(weight_percentages,digits=1), nsmall=1)\nwp = weight_percentages\/100\n#Forest plot box sizes on weihts\nwbox=c(NA,NA,(k\/4)*sqrt(wp)\/sum(sqrt(wp)),0.5,0,NA)\t\n#wbox=c(NA,NA,(k\/5)*sqrt(wp)\/sum(sqrt(wp)),0.5,0,NA)\n###### tau-squared\n# RD \u306e tau^2\ntau_squared_RD_samples <- posterior_samples$tau_squared_RD\n#mean_tau_squared_RD <- mean(tau_squared_RD_samples)\n#ci_tau_squared_RD <- quantile(tau_squared_RD_samples, probs = c(0.025, 0.975))\n#mode and HDI (High Density Interval) logOR tau-squared\ndens=density(tau_squared_RD_samples)\nmean_tau_squared_RD=dens$x[which.max(dens$y)]\t#mode\nci_tau_squared_RD=hdi(tau_squared_RD_samples)\t#High Density Interval (HDI) lower and upper.\n\nprint(paste(\"Estimated tau-squared (RD, mode):\", mean_tau_squared_RD))\nprint(paste(\"95% High Density Interval for tau-squared (RD):\", ci_tau_squared_RD[1], \"-\", ci_tau_squared_RD[2]))\n\n#Forest plot by forestplot\n#setting fs for cex\nfs=1\nif(k>20){fs=round((1-0.02*(k-20)),digits=1)}\n\nm=c(NA,NA,rd,risk_difference,risk_difference,NA)\nlw=c(NA,NA,rd_lw,conf_int_rd[1],pi_rd[1],NA)\nup=c(NA,NA,rd_up,conf_int_rd[2],pi_rd[2],NA)\n\nhete1=\"\"\nhete2=paste(\"I2=\",round(mean_I_squared_RD, 2),\"%\",sep=\"\")\n\nhete3=paste(\"tau2=\",format(round(mean_tau_squared_RD,digits=4),nsmall=4),sep=\"\")\n\n#hete2=paste(\"I2=\",round(mean_I_squared_RD, 2), \"%\\n\",\" (\",round(ci_I_squared_RD[1], 2),\" ~ \",\n#round(ci_I_squared_RD[2], 2),\")\",sep=\"\")\n\n#hete3=paste(\"tau2=\",format(round(mean_tau_squared_RD,digits=4),nsmall=4), \"\\n\",\n#\" (\",format(round(ci_tau_squared_RD[1],digits=4),nsmall=4),\" ~ \",\n#format(round(ci_tau_squared_RD[2],digits=4),nsmall=4),\")  \",sep=\"\")\n\nhete4=\"\"\nhete5=paste(\"p=\",p_val_rd, sep=\"\")\nhete6=paste(label,prob_direct_rd,sep=\"\")\n\nau=study\nsl=c(NA,toString(labe_study),as.vector(au),\"Summary Estimate\",\"Prediction Interval\",NA)\n\nncl=c(labe_cont,\"Number\",nc,NA,NA,hete1)\nrcl=c(NA,labe_outc,rc,NA,NA,hete2)\n\nntl=c(labe_int,\"Number\",nt,NA,NA,hete3)\nrtl=c(NA,labe_outc,rt,NA,NA,hete4)\n\nspac=c(\"    \",NA,rep(NA,k),NA,NA,NA)\nml=c(NA,labe_em,format(round(rd,digits=3),nsmall=3),format(round(risk_difference,digits=3),nsmall=3),NA,hete5)\nll=c(NA,\"95% CI lower\",format(round(rd_lw,digits=3),nsmall=3),format(round(conf_int_rd[1],digits=3),nsmall=3),format(round(pi_rd[1],digits=3),nsmall=3),hete6)\nul=c(NA,\"95% CI upper\",format(round(rd_up,digits=3),nsmall=3),format(round(conf_int_rd[2],digits=3),nsmall=3),format(round(pi_rd[2],digits=3),nsmall=3),NA)\nwpcl=c(NA,\"Weight(%)\",wpc,100,NA,NA)\nll=as.vector(ll)\nul=as.vector(ul)\nsum=c(TRUE,TRUE,rep(FALSE,k),TRUE,FALSE,FALSE)\nzerov=0\nlabeltext=cbind(sl,ntl,rtl,ncl,rcl,spac,ml,ll,ul,wpcl)\nhlines=list(\"3\"=gpar(lwd=1,columns=1:11,col=\"grey\"))\ndev.new()\nplot(forestplot(labeltext,mean=m,lower=lw,upper=up,is.summary=sum,graph.pos=7,\nzero=zerov,hrzl_lines=hlines,xlab=toString(labe_em),txt_gp=fpTxtGp(ticks=gpar(cex=fs),\nxlab=gpar(cex=fs),cex=fs),xticks.digits=2,vertices=TRUE,graphwidth=unit(50,\"mm\"),colgap=unit(3,\"mm\"),\nboxsize=wbox,\nlineheight=\"auto\",xlog=FALSE,new_page=FALSE))\n}\n####\n###Pint indiv estimates with CI and copy to clipboard.\nnk=k+2\nefs=rep(\"\",nk)\nefeci=rep(\"\",nk)\nfor(i in 1:nk){\nefs[i]=ml[i+2]\nefeci[i]=paste(ll[i+2],\"~\",ul[i+2])\n}\nefestip=data.frame(cbind(efs,efeci))\nefestip[nk,1]=\"\"\nprint(efestip)\nwrite.table(efestip,\"clipboard\",sep=\"\\t\",row.names=FALSE,col.names=FALSE)\nprint(\"The estimates of each study and the summary estimate are in the clipboard.\")\n#############################################<\/code><\/pre>\n<\/div><\/details>\n\n\n\n
\"\"<\/figure>\n\n\n\n
\u51fa\u529b\u3055\u308c\u305fForest plot\u3067\u3059\u304c\u3001\u78ba\u4fe1\u533a\u9593\u3060\u3051\u3067\u306a\u304f\u4e88\u6e2c\u533a\u9593\u3082\u51fa\u529b\u3057\u307e\u3059\u3002Forest plot\u304c\u51fa\u529b\u3055\u308c\u305f\u6642\u70b9\u3067\u5404\u7814\u7a76\u306e\u52b9\u679c\u63a8\u5b9a\u5024\u300195\uff05\u78ba\u4fe1\u533a\u9593\u3001\u7d71\u5408\u5024\u3068\u305d\u306e95\uff05\u78ba\u4fe1\u533a\u9593\u3001\u4e88\u6e2c\u533a\u9593\u306e\u5024\u3092\u30b3\u30f3\u30bd\u30fc\u30eb\u306b\u51fa\u529b\u3059\u308b\u3068\u3068\u3082\u306b\u3001\u30af\u30ea\u30c3\u30d7\u30dc\u30fc\u30c9\u306b\u30b3\u30d4\u30fc\u3057\u307e\u3059\u306e\u3067\u3001Excel\u306e\u8a55\u4fa1\u30b7\u30fc\u30c8\u306a\u3069\u306b\u7d50\u679c\u306e\u5024\u3092\u8cbc\u308a\u4ed8\u3051\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u3002<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
\u5fc5\u8981\u306a\u30d1\u30c3\u30b1\u30fc\u30b8\u304c\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u3055\u308c\u3066\u3044\u306a\u3044\u5834\u5408\u3001\u4ee5\u4e0b\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u5148\u306b\u5b9f\u884c\u3057\u3066\u304f\u3060\u3055\u3044\u3002\u672a\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u306e\u30d1\u30c3\u30b1\u30fc\u30b8\u304c\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u3055\u308c\u307e\u3059\u3002metafor\u3082\u542b\u3081\u307e\u3057\u305f\u3002\u307e\u305f\u3001R\u3067Stan\u3092\u52d5\u304b\u3059\u305f\u3081\u306b\u306f\u3001RTools\u306e\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u3082\u5fc5\u8981\u306b\u306a\u308a\u307e\u3059\u3002\u3053\u3061\u3089\u304b\u3089\u3001R\u306e\u30d0\u30fc\u30b8\u30e7\u30f3\u306b\u5408\u308f\u305b\u305f\u30d0\u30fc\u30b8\u30e7\u30f3\u306eRTools\u3092\u30c0\u30a6\u30f3\u30ed\u30fc\u30c9\u3057\u3066\u30a4\u30f3\u30b9\u30c8\u30fc\u30eb\u3057\u3066\u304f\u3060\u3055\u3044\uff08430MB \u304f\u3089\u3044\u3042\u308b\u5927\u304d\u306a\u30d7\u30ed\u30b0\u30e9\u30e0\u3067\u3059\uff09\u3002URL https:\/\/cran.r-project.org\/bin\/windows\/Rtools\/<\/a> \u8a73\u7d30\u306f\u524d\u306e\u6295\u7a3f<\/a>\u3092\u53c2\u7167\u3057\u3066\u304f\u3060\u3055\u3044\u3002<\/p>\n\n\n\n
####Install packages for the first time####\npackneed=c(\"rstan\",\"Rcpp\",\"ggplot2\",\"coda\",\"dplyr\",\"forestplot\",\"bayesplot\",\"HDInterval\",\"metafor\");current=installed.packages();addpack=setdiff(packneed,rownames(current));url=\"https:\/\/ftp.yz.yamagata-u.ac.jp\/pub\/cran\/\";if(length(addpack)>0){install.packages(addpack,repos=url)};if(length(addpack)==0){print(\"Already installed.\")}<\/code><\/pre>\n\n\n\n
Funnel plot\u4f5c\u6210\u306e\u305f\u3081\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u4f5c\u6210\u3057\u307e\u3057\u305f\u3002Forest plot\u304c\u51fa\u529b\u3055\u308c\u305f\u5f8c\u306b\u5b9f\u884c\u3059\u308b\u3068Funnel plot\u304c\u51fa\u529b\u3055\u308c\u307e\u3059\u3002\u307e\u305f\u3001\u7d71\u8a08\u5b66\u7684\u306a\u975e\u5bfe\u79f0\u6027\u306e\u691c\u5b9a\u6cd5\u3067\u3042\u308b\u3001Begg\u306e\u691c\u5b9a\u3001Egger\u306e\u691c\u5b9a\u3082\u5b9f\u884c\u3057\u307e\u3059\u3002\u6700\u5f8c\u307e\u3067\u5b9f\u884c\u3059\u308b\u3068Funnel plot\u5185\u306b\u3053\u308c\u3089\u691c\u5b9a\u306e\u7d50\u679c\u3082\u66f8\u304d\u51fa\u3057\u307e\u3059\u3002Funnel plot\u3060\u3051\u3067\u5341\u5206\u306a\u5834\u5408\u306f\u3001\u5f8c\u534a\u306eAsymmetry test with Egger and Begg’s tests.\u306e\u884c\u4ee5\u4e0b\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u306f\u5b9f\u884c\u3059\u308b\u5fc5\u8981\u304c\u3042\u308a\u307e\u305b\u3093\u3002\u307e\u305f\u3001Funnel plot\u306f\u7814\u7a76\u6570\u304c10\uff08\u5c11\u306a\u304f\u3068\u30825\uff09\u4ef6\u4ee5\u4e0a\u306e\u5834\u5408\u306b\u610f\u5473\u304c\u3042\u308b\u3068\u3055\u308c\u3066\u3044\u307e\u3059\u304c\u3001\u3053\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3067\u306f\u6570\u306b\u3088\u308b\u5236\u9650\u306f\u3057\u3066\u3044\u307e\u305b\u3093\u306e\u3002<\/p>\n\n\n\n
Funnel plot\u4f5c\u6210\u306e\u305f\u3081\u306e\u30b9\u30af\u30ea\u30d7\u30c8<\/summary>
\n
########################\nlibrary(metafor)\n\n#Plot funnel plot.\nif(em == \"RR\"){\nyi = log(rr)\nvi=1\/weights\nsei = sqrt(1\/weights)\nmu = log(risk_ratio)\ndev.new(width=7,height=7)\nfunnel(yi, sei = sei, level = c(95), refline = mu,xlab = \"log(RR)\", ylab = \"Standard Error\")\n}\n\nif(em == \"OR\"){\nyi = log(or)\nvi=1\/weights\nsei = sqrt(1\/weights)\nmu = log(odds_ratio)\ndev.new(width=7,height=7)\nfunnel(yi, sei = sei, level = c(95), refline = mu,xlab = \"log(OR)\", ylab = \"Standard Error\")\n}\n\nif(em == \"RD\"){\nyi = rd\nvi=1\/weights\nsei = sqrt(1\/weights)\nmu = risk_difference\ndev.new(width=7,height=7)\nfunnel(yi, sei = sei, level = c(95), refline = mu,xlab = \"RD\", ylab = \"Standard Error\")\n}\n\n#Asymmetry test with Egger and Begg's tests.\negger=regtest(yi, sei=sei,model=\"lm\", ret.fit=FALSE)\nbegg=ranktest(yi, vi)\n\n#Print the results to the console.\nprint(\"Egger's test:\")\nprint(egger)\nprint(\"Begg's test\")\nprint(begg)\n\n#Add Begg and Egger to the plot.\nfsfn=1\nem=toString(exdat$label[6])\noutyes=toString(exdat$label[4])\nfunmax=par(\"usr\")[3]-par(\"usr\")[4]\ngyou=funmax\/12\nfxmax=par(\"usr\")[2]-(par(\"usr\")[2]-par(\"usr\")[1])\/40\ntext(fxmax,gyou*0.5,\"Begg's test\",pos=2,cex=fsfn)\nkentau=toString(round(begg$tau,digits=3))\ntext(fxmax,gyou*1.2,paste(\"Kendall's tau=\",kentau,sep=\"\"),pos=2,cex=fsfn)\nkenp=toString(round(begg$pval,digits=5))\ntext(fxmax,gyou*1.9,paste(\"p=\",kenp,sep=\"\"),pos=2,cex=fsfn)\n\ntext(fxmax,gyou*3.5,\"Egger's test\",pos=2,cex=fsfn)\ntstat=toString(round(egger$zval,digits=3))\ntext(fxmax,gyou*4.2,paste(\"t statistic=\",tstat,sep=\"\"),pos=2,cex=fsfn)\ntstap=toString(round(egger$pval,digits=5))\ntext(fxmax,gyou*5.1,paste(\"p=\",tstap,sep=\"\"),pos=2,cex=fsfn)\n#####################################<\/code><\/pre>\n<\/div><\/details>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
\u4ee5\u524d\u306e\u6295\u7a3f\u300cStan\u3092\u4f7f\u3046\u30d9\u30a4\u30b8\u30a2\u30f3\u30a2\u30d7\u30ed\u30fc\u30c1\u306e\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\uff1aCopilot\u306b\u805e\u304d\u306a\u304c\u3089\u300d\u3067R\u304b\u3089rstan\u3092\u4ecb\u3057\u3066Stan\u3092\u52d5\u304b\u3057\u3001\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3092\u5b9f\u884c\u3059\u308b\u307e\u3067\u306e\u904e\u7a0b\u3092\u7d39\u4ecb\u3057\u307e\u3057\u305f\u3002 \u4eca\u56de\u3001\u30ea\u30b9\u30af\u6bd4\u3001\u30aa\u30c3\u30ba\u6bd4\u3001\u30ea\u30b9\u30af … “Stan\u3092\u4f7f\u3046\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u306e\u30b3\u30fc\u30c9\uff1a\u4e8c\u5024\u5909\u6570\u30a2\u30a6\u30c8\u30ab\u30e0RR, OR, RD\u7528” \u306e<\/span>\u7d9a\u304d\u3092\u8aad\u3080<\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[32,2],"tags":[],"class_list":["post-1885","post","type-post","status-publish","format-standard","hentry","category-statistics","category-sr"],"_links":{"self":[{"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/posts\/1885","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=1885"}],"version-history":[{"count":7,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/posts\/1885\/revisions"}],"predecessor-version":[{"id":1906,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/posts\/1885\/revisions\/1906"}],"wp:attachment":[{"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1885"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1885"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1885"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}