{"id":1893,"date":"2025-10-30T22:15:01","date_gmt":"2025-10-30T13:15:01","guid":{"rendered":"https:\/\/info.zanet.biz\/?p=1893"},"modified":"2025-11-01T16:04:09","modified_gmt":"2025-11-01T07:04:09","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-2","status":"publish","type":"post","link":"https:\/\/info.zanet.biz\/?p=1893","title":{"rendered":"Stan\u3092\u4f7f\u3046\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u306e\u30b3\u30fc\u30c9\uff1a\u9023\u7d9a\u5909\u6570\u30a2\u30a6\u30c8\u30ab\u30e0MD\u3001SMD\u7528"},"content":{"rendered":"\n

\u9023\u7d9a\u5909\u6570\u30a2\u30a6\u30c8\u30ab\u30e0\u3067\u5e73\u5747\u5024\u5deeMean Difference, MD\u307e\u305f\u306f\u6a19\u6e96\u5316\u5e73\u5747\u5024\u5deeStandardized Mean Difference (SMD)\u306e\u7d71\u5408\u5024\u306895\uff05\u78ba\u4fe1\u533a\u9593Credible Interval\u3092\u3001Stan\u3092\u7528\u3044\u3066\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u3067\u6c42\u3081\u3001Forest plot\u3092\u4f5c\u6210\u3059\u308b\u30b3\u30fc\u30c9\u3067\u3059\u3002<\/p>\n\n\n\n

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

\u30c7\u30fc\u30bf\u306fExcel\u3067\u56f3\u306e\u3088\u3046\u306b\u6e96\u5099\u3057\u307e\u3059\u3002label\u306f\u7814\u7a76\u540d\u306e\u30e9\u30d9\u30eb\u3001\u5bfe\u7167\u3001\u4ecb\u5165\u306e\u30e9\u30d9\u30eb\u3001\u30a2\u30a6\u30c8\u30ab\u30e0\u3001\u52b9\u679c\u6307\u6a19\u306e\u30bf\u30a4\u30d7\u3001\u305d\u306e\u7565\u79f0\u3067\u3059\u3002\u30bb\u30ebI9\u304cMD\u306e\u5834\u5408\u306f\u3001\u5e73\u5747\u5024\u5dee\u3001SMD\u306e\u5834\u5408\u306f\u6a19\u6e96\u5316\u5e73\u5747\u5024\u5dee\u3092\u8a08\u7b97\u3057\u307e\u3059\u3002\u3053\u306e\u4f8b\u3067\u306fSMD\u3092\u8a08\u7b97\u3057\u307e\u3059\u3002<\/p>\n\n\n\n

R\u3092\u8d77\u52d5\u3057\u3066\u304a\u304d\u3001\u30a8\u30c7\u30a3\u30bf\u30fc\u306b\u30b3\u30d4\u30fc\u3057\u305f\u4e0b\u8a18\u30b3\u30fc\u30c9\u3092\u8cbc\u308a\u4ed8\u3051\u3066\u3001\u6700\u521d\u306e5\u3064\u306e\u30e9\u30a4\u30d6\u30e9\u30ea\u306e\u8aad\u307f\u8fbc\u307f\u3092\u5b9f\u884c\u3057\u3066\u304a\u304d\u307e\u3059\u3002Excel\u306b\u623b\u3063\u3066\u3001\u30bb\u30ebB3\u304b\u3089I10\u307e\u3067\u306e\u7bc4\u56f2\u3092\u9078\u629e\u3057\u3066\u3001\u30b3\u30d4\u30fc\u64cd\u4f5c\uff08Ctrl+C)\u3092\u884c\u3063\u3066\u3001R\u306b\u623b\u308a\u3001exdato = read.delim(“clipboard”,sep=”\\t”,header=TRUE)\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u5b9f\u884c\u3057\u307e\u3059\u3002\u7d9a\u3044\u3066\u3001\u305d\u308c\u4ee5\u4e0b\u306e\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u5168\u90e8\u5b9f\u884c\u3055\u305b\u308b\u3068\u3001Forest plot\u304c\u51fa\u529b\u3055\u308c\u307e\u3059\u3002\u52b9\u679c\u6307\u6a19\u306e\u70b9\u63a8\u5b9a\u5024\u306895\uff05\u78ba\u4fe1\u533a\u9593\u3001\u304a\u3088\u3073\u4e88\u6e2c\u533a\u9593\u3092\u63d0\u793a\u3057\u307e\u3059\u3002<\/p>\n\n\n\n

Markov Chain Monte Carlo (MCMC)\u30b7\u30df\u30e5\u30ec\u30fc\u30b7\u30e7\u30f3\u3092\u884c\u3046\u306e\u3067\u3001\u5c11\u3057\u6642\u9593\u304c\u304b\u304b\u308a\u307e\u3059\u3002chains = 4, warmup = 10000, iter = 20000\u306b\u8a2d\u5b9a\u3057\u3066\u3044\u308b\u306e\u3067\u3001\u8a0840000\u500b\u304c\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u3055\u308c\u307e\u3059\u3002<\/p>\n\n\n\n

Forest plot\u304c\u51fa\u529b\u3055\u308c\u305f\u6642\u70b9\u3067\u3001\u5404\u7814\u7a76\u306e\u52b9\u679c\u63a8\u5b9a\u5024\u300195\uff05\u78ba\u4fe1\u533a\u9593\u3001\u7d71\u5408\u5024\u306895\uff05\u78ba\u4fe1\u533a\u9593\u3001\u4e88\u6e2c\u533a\u9593\u306e\u5024\u3092\u30af\u30ea\u30c3\u30d7\u30dc\u30fc\u30c9\u306b\u30b3\u30d4\u30fc\u3057\u3066\u3044\u307e\u3059\u306e\u3067\u3001Excel\u8a55\u4fa1\u30b7\u30fc\u30c8\u306a\u3069\u306b\u8cbc\u308a\u4ed8\u3051\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u3002<\/p>\n\n\n\n

\u5e73\u5747\u5024\u5deeMD\u3001\u6a19\u6e96\u5316\u5e73\u5747\u5024\u5deeSMD\u306e\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u30b3\u30fc\u30c9<\/summary>
\n
#####MD and SMD Bayesian meta-analysis with Stan, rstan######\n# R code for data preparation\nlibrary(rstan)\nlibrary(bayesplot)\nlibrary(dplyr)\nlibrary(HDInterval)\nlibrary(forestplot)\nlibrary(metafor)\n\n# Get data via clipboard from Excel sheet.\nexdato = read.delim(\"clipboard\",sep=\"\\t\",header=TRUE)\nlabem = exdato[,\"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\nexdat=na.omit(exdato)\n# Number of studies\nK <- nrow(exdat)\n\n# Calculate individual study effect sizes (MD and SMD) and their variances\n# For Mean Difference (MD)\nexdat$yi_md <- exdat$m1i - exdat$m2i\nexdat$vi_md <- (exdat$sd1i^2 \/ exdat$n1i) + (exdat$sd2i^2 \/ exdat$n2i)\nexdat$sei_md <- sqrt(exdat$vi_md)\n\n# For Standardized Mean Difference (SMD) - using Hedges' g approximation\n# Calculate pooled standard deviation\nexdat$sd_pooled <- with(exdat, sqrt(((n1i - 1) * sd1i^2 + (n2i - 1) * sd2i^2) \/ (n1i + n2i - 2)))\nexdat$yi_smd <- (exdat$m1i - exdat$m2i) \/ exdat$sd_pooled\n\n# Calculate variance for SMD (Hedges' g)\nexdat$vi_smd <- with(exdat, ((n1i + n2i) \/ (n1i * n2i)) + (yi_smd^2 \/ (2 * (n1i + n2i))))\nexdat$sei_smd <- sqrt(exdat$vi_smd)\n\nif(em==\"MD\"){\n# Stan data for MD\nstan_data_md <- list(\n  K = K,\n  yi = exdat$yi_md,\n  sei = exdat$sei_md\n)\n}\nif(em==\"SMD\"){\n# Stan data for SMD\nstan_data_smd <- list(\n  K = K,\n  yi = exdat$yi_smd,\n  sei = exdat$sei_smd\n)\n}\n######\nif(em==\"MD\"){\n## Stan code for Mean Difference (MD) random-effects model\nstan_md_code <- \"\ndata {\n  int<lower=0> K;        \/\/ number of studies\n  array[K] real yi;      \/\/ observed effect sizes\n  array[K] real<lower=0> sei; \/\/ standard errors of effect sizes\n}\n\nparameters {\n  real mu;               \/\/ overall mean effect\n  real<lower=0> tau;     \/\/ between-study standard deviation\n  array[K] real delta;   \/\/ true effect size for each study\n}\n\ntransformed parameters {\n  real<lower=0> tau_squared; \/\/ between-study variance\n  real<lower=0,upper=1> I_squared; \/\/ I-squared statistic\n\n  tau_squared = square(tau);\n\n  \/\/ Calculate I-squared\n  \/\/ Average sampling variance (approximation)\n  real var_sampling_avg = mean(square(sei));\n  I_squared = tau_squared \/ (tau_squared + var_sampling_avg);\n}\n\nmodel {\n  \/\/ Priors\n  mu ~ normal(0, 10);      \/\/ weakly informative prior for overall mean\n  tau ~ normal(0, 1);      \/\/ weakly informative prior for tau (sd of random effects)\n\n  \/\/ Likelihood\n  delta ~ normal(mu, tau); \/\/ true effect sizes drawn from a normal distribution\n  yi ~ normal(delta, sei); \/\/ observed effect sizes are normally distributed around true effect sizes\n}\n\ngenerated quantities {\n  real prediction_interval_lower;\n  real prediction_interval_upper;\n  real new_delta; \/\/ New true effect size from the meta-analytic distribution\n\n  \/\/ Prediction interval for a new study's true effect size\n  new_delta = normal_rng(mu, tau);\n  \n  \/\/ To get the interval, we sample many 'new_delta' and take quantiles later\n  \/\/ For simplicity, we can also directly calculate it if mu and tau are well-behaved.\n  \/\/ However, for a true Bayesian prediction interval, it's best to sample.\n  \/\/ Here, for outputting directly from Stan, we can use quantiles of posterior samples of mu and tau\n  \/\/ and combine them, or draw samples as below.\n  \/\/ For a direct 95% interval, we can compute it from mu and tau:\n  prediction_interval_lower = mu - 1.96 * tau;\n  prediction_interval_upper = mu + 1.96 * tau;\n}\n\"\n}\n#####\nif(em == \"SMD\"){\n## Stan code for Standardized Mean Difference (SMD) random-effects model\nstan_smd_code <- \"\ndata {\n  int<lower=0> K;        \/\/ number of studies\n  array[K] real yi;      \/\/ observed effect sizes (SMD)\n  array[K] real<lower=0> sei; \/\/ standard errors of effect sizes\n}\n\nparameters {\n  real mu;               \/\/ overall mean effect (SMD)\n  real<lower=0> tau;     \/\/ between-study standard deviation\n  array[K] real delta;   \/\/ true effect size for each study\n}\n\ntransformed parameters {\n  real<lower=0> tau_squared; \/\/ between-study variance\n  real<lower=0,upper=1> I_squared; \/\/ I-squared statistic\n\n  tau_squared = square(tau);\n\n  \/\/ Calculate I-squared\n  \/\/ Average sampling variance (approximation)\n  real var_sampling_avg = mean(square(sei));\n  I_squared = tau_squared \/ (tau_squared + var_sampling_avg);\n}\n\nmodel {\n  \/\/ Priors\n  mu ~ normal(0, 1);      \/\/ weakly informative prior for overall mean (SMD is typically smaller scale)\n  tau ~ normal(0, 0.5);   \/\/ weakly informative prior for tau (sd of random effects, typically smaller for SMD)\n\n  \/\/ Likelihood\n  delta ~ normal(mu, tau); \/\/ true effect sizes drawn from a normal distribution\n  yi ~ normal(delta, sei); \/\/ observed effect sizes are normally distributed around true effect sizes\n}\n\ngenerated quantities {\n  real prediction_interval_lower;\n  real prediction_interval_upper;\n  real new_delta; \/\/ New true effect size from the meta-analytic distribution\n\n  \/\/ Prediction interval for a new study's true effect size\n  new_delta = normal_rng(mu, tau);\n  prediction_interval_lower = mu - 1.96 * tau;\n  prediction_interval_upper = mu + 1.96 * tau;\n}\n\"\n}\n#####\n# R code to run Stan models and extract results\nif(em==\"MD\"){\n# Compile the MD Stan model\nstan_md_model <- stan_model(model_code = stan_md_code)\n\n# Run the MD model\nfit_md <- sampling(stan_md_model,\n                   data = stan_data_md,\n                   chains = 4,             # number of MCMC chains\n                   iter = 20000,            # number of iterations per chain\n                   warmup = 10000,          # warmup iterations\n                   thin = 1,               # thinning rate\n                   seed = 1234,            # for reproducibility\n                   control = list(adapt_delta = 0.95, max_treedepth = 15)) # improve sampling\n\nprint(fit_md, pars = c(\"mu\", \"tau\", \"tau_squared\", \"I_squared\", \"prediction_interval_lower\", \"prediction_interval_upper\"),\n      probs = c(0.025, 0.5, 0.975))\n\n# Plotting (optional)\n# dev.new();mcmc_dens(fit_md, pars = c(\"mu\", \"tau\", \"tau_squared\", \"I_squared\"))\n# dev.new();mcmc_trace(fit_md, pars = c(\"mu\", \"tau\"))\n}\n###\nif(em==\"SMD\"){\n# Compile the SMD Stan model\nstan_smd_model <- stan_model(model_code = stan_smd_code)\n\n# Run the SMD model\nfit_smd <- sampling(stan_smd_model,\n                    data = stan_data_smd,\n                    chains = 4,\n                    iter = 20000,\n                    warmup = 10000,\n                    thin = 1,\n                    seed = 1234,\n                    control = list(adapt_delta = 0.95, max_treedepth = 15))\n\nprint(fit_smd, pars = c(\"mu\", \"tau\", \"tau_squared\", \"I_squared\", \"prediction_interval_lower\", \"prediction_interval_upper\"),\n      probs = c(0.025, 0.5, 0.975))\n\n# Plotting (optional)\n# dev.new();mcmc_dens(fit_smd, pars = c(\"mu\", \"tau\", \"tau_squared\", \"I_squared\"))\n# dev.new();mcmc_trace(fit_smd, pars = c(\"mu\", \"tau\"))\n}\n#####\nif(em==\"MD\"){\n#Posterior samples of MD\nposterior_samples=extract(fit_md)\n#md meand difference\nmu_md=mean(posterior_samples$mu)\nmu_md_ci=quantile(posterior_samples$mu,probs=c(0.025,0.975))\nmu_md_lw=mu_md_ci[1]\nmu_md_up=mu_md_ci[2]\np_val_md = round(2*min(mean(posterior_samples$mu >0), mean(posterior_samples$mu <0)), 6)\t#Bayesian p-value for summary MD\nif(mu_md>0){\nprob_direct_md = round(mean(posterior_samples$mu > 0),6)\nlabel = \"p(MD>0)=\"\n}else{\nprob_direct_md = round(mean(posterior_samples$mu < 0),6)\nlabel = \"p(MD<0)=\"\n}\n#md tau-squared\n#if normal distribution.\n#mu_md_tau_squared=mean(posterior_samples$tau_squared)\n#mu_md_tau_squared_ci=quantile(posterior_samples$tau_squared,probs=c(0.025,0.975))\n#mu_md_tau_squared_lw=mu_md_tau_squared_ci[1]\n#mu_md_tau_squared_up=mu_md_tau_squared_ci[2]\n#if skewed distributions.\ndens=density(posterior_samples$tau_squared)\nmu_md_tau_squared=dens$x[which.max(dens$y)]\t#mode\nmu_md_tau_squared_lw=hdi(posterior_samples$tau_squared)[1]\t#High Density Interval (HDI) lower\nmu_md_tau_squared_up=hdi(posterior_samples$tau_squared)[2]\t#High Density Interval (HDI) upper\n\n#md I-squared\n#if normal distribution.\n#mu_md_I_squared=mean(posterior_samples$I_squared)\n#mu_md_I_squared_ci=quantile(posterior_samples$I_squared,probs=c(0.025,0.975))\n#mu_md_I_squared_lw=mu_md_I_squared_ci[1]\n#mu_md_I_squared_up=mu_md_I_squared_ci[2]\n#if skewed distributions.\ndens=density(posterior_samples$I_squared)\nmu_md_I_squared=dens$x[which.max(dens$y)]\t#mode\nmu_md_I_squared_lw=hdi(posterior_samples$I_squared)[1]\t#High Density Interval (HDI) lower\nmu_md_I_squared_up=hdi(posterior_samples$I_squared)[2]\t#High Density Interval (HDI) upper\n\n#md Prediction Interval\nmu_new_md=mean(posterior_samples$new_delta)\nmu_new_md_ci=quantile(posterior_samples$new_delta,probs=c(0.025,0.975))\nmu_md_pi_lw=median(posterior_samples$prediction_interval_lower)\nmu_md_pi_up=median(posterior_samples$prediction_interval_upper)\n}\n####\nif(em==\"SMD\"){\n#Posterior samples of SMD\nposterior_samples=extract(fit_smd)\n#smd meand difference\nmu_smd=mean(posterior_samples$mu)\nmu_smd_ci=quantile(posterior_samples$mu,probs=c(0.025,0.975))\nmu_smd_lw=mu_smd_ci[1]\nmu_smd_up=mu_smd_ci[2]\np_val_smd = round(2*min(mean(posterior_samples$mu >0), mean(posterior_samples$mu <0)), 6)\t#Bayesian p-value for summary MD\nif(mu_smd>0){\nprob_direct_smd = round(mean(posterior_samples$mu > 0),6)\nlabel = \"p(SMD>0)=\"\n}else{\nprob_direct_smd = round(mean(posterior_samples$mu < 0),6)\nlabel = \"p(SMD<0)=\"\n}\n#smd tau-squared\n#if normal distribution.\n#mu_smd_tau_squared=mean(posterior_samples$tau_squared)\t#mean\n#mu_smd_tau_squared_ci=quantile(posterior_samples$tau_squared,probs=c(0.025,0.975))\n#mu_smd_tau_squared_lw=mu_smd_tau_squared_ci[1]\n#mu_smd_tau_squared_up=mu_smd_tau_squared_ci[2]\n#if skewed distributions.\ndens=density(posterior_samples$tau_squared)\nmu_smd_tau_squared=dens$x[which.max(dens$y)]\t#mode\nmu_smd_tau_squared_lw=hdi(posterior_samples$tau_squared)[1]\t#High Density Interval (HDI) lower\nmu_smd_tau_squared_up=hdi(posterior_samples$tau_squared)[2]\t#High Density Interval (HDI) upper\n\n#smd I-squared\n#if normal distribution.\n#mu_smd_I_squared=mean(posterior_samples$I_squared)\n#mu_smd_I_squared_ci=quantile(posterior_samples$I_squared,probs=c(0.025,0.975))\n#mu_smd_I_squared_lw=mu_smd_I_squared_ci[1]\n#mu_smd_I_squared_up=mu_smd_I_squared_ci[2]\n#if skewed distributions.\ndens=density(posterior_samples$I_squared)\nmu_smd_I_squared=dens$x[which.max(dens$y)]\t#mode\nmu_smd_I_squared_lw=hdi(posterior_samples$I_squared)[1]\t#High Density Interval (HDI) lower\nmu_smd_I_squared_up=hdi(posterior_samples$I_squared)[2]\t#High Density Interval (HDI) upper\n\n#smd Prediction Interval\nmu_new_smd=mean(posterior_samples$new_delta)\nmu_new_smd_ci=quantile(posterior_samples$new_delta,probs=c(0.025,0.975))\nmu_smd_pi_lw=median(posterior_samples$prediction_interval_lower)\nmu_smd_pi_up=median(posterior_samples$prediction_interval_upper)\n}\n#####\nif(em==\"MD\"){\n#\u5404\u7814\u7a76\u306emd\u306895%CI\nmd=rep(0,K)\nmd_lw=md\nmd_up=md\nfor(i in 1:K){\nmd[i]=mean(posterior_samples$delta[,i])\nmd_lw[i]=quantile(posterior_samples$delta[,i],probs=0.025)\nmd_up[i]=quantile(posterior_samples$delta[,i],probs=0.975)\n}\n}\nif(em==\"SMD\"){\n#\u5404\u7814\u7a76\u306esmd\u306895%CI\nsmd=rep(0,K)\nsmd_lw=smd\nsmd_up=smd\nfor(i in 1:K){\nsmd[i]=mean(posterior_samples$delta[,i])\nsmd_lw[i]=quantile(posterior_samples$delta[,i],probs=0.025)\nsmd_up[i]=quantile(posterior_samples$delta[,i],probs=0.975)\n}\n}\n#########\n#Forest Plot\nif(em==\"MD\"){\n#MD\n# \u91cd\u307f\u3092\u683c\u7d0d\u3059\u308b\u30d9\u30af\u30c8\u30eb\u3092\u521d\u671f\u5316\nweights <- numeric(K)\nweight_percentages <- numeric(K)\n\nfor (i in 1:K) {\n  # i\u756a\u76ee\u306e\u7814\u7a76\u306emd\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u5024\u3092\u53d6\u5f97\n  md_i_samples <- posterior_samples$delta[, 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_md_i <- var(md_i_samples)\n  # \u91cd\u307f\u306f\u5206\u6563\u306e\u9006\u6570\n  weights[i] <- 1 \/ variance_of_md_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:K,\n  Effective_Weight = weights,\n  Effective_Weight_Percentage = weight_percentages\n)\nprint(results_weights_posterior_var)\nk=K\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\n#setting fs for cex\nfs=1\nif(k>20){fs=round((1-0.02*(k-20)),digits=1)}\n\nm=c(NA,NA,md,mu_md,mu_new_md,NA)\nlw=c(NA,NA,md_lw,mu_md_lw,mu_md_pi_lw,NA)\nup=c(NA,NA,md_up,mu_md_up,mu_md_pi_up,NA)\n\nhete1=\"\"\nhete2=paste(\"I2=\",round(100*mu_md_I_squared, 2),\"%\",sep=\"\")\n\nhete3=paste(\"tau2=\",format(round(mu_md_tau_squared,digits=4),nsmall=4),sep=\"\")\nhete4=\"\"\nhete5=paste(\"p=\",p_val_md, sep=\"\")\nhete6=paste(label,prob_direct_md, sep=\"\")\n\nrck=rep(NA,K)\nrtk=rck\nfor(i in 1:K)\n{\nrtk[i]=paste(exdat$m1i[i],\" (\",exdat$sd1i[i],\")\",sep=\"\")\nrck[i]=paste(exdat$m2i[i],\" (\",exdat$sd2i[i],\")\",sep=\"\")\n}\n\nau=exdat$author\nsl=c(NA,toString(labe_study),as.vector(au),\"Summary Estimate\",\"Prediction Interval\",NA)\nnc=exdat$n2i\nncl=c(labe_cont,\"Number\",nc,NA,NA,hete1)\nrcl=c(labe_outc,\"Mean (SD)\",rck,NA,NA,hete2)\n\nnt=exdat$n1i\nntl=c(labe_int,\"Number\",nt,NA,NA,hete3)\nrtl=c(labe_outc,\"Mean (SD)\",rtk,NA,NA,hete4)\n\nspac=c(\"    \",NA,rep(NA,k),NA,NA,NA)\nml=c(NA,labe_em,format(round(md,digits=3),nsmall=3),format(round(mu_md,digits=3),nsmall=3),NA,hete5)\nll=c(NA,\"95% CI lower\",format(round(md_lw,digits=3),nsmall=3),format(round(mu_md_lw,digits=3),nsmall=3),format(round(mu_md_pi_lw,digits=3),nsmall=3),hete6)\nul=c(NA,\"95% CI upper\",format(round(md_up,digits=3),nsmall=3),format(round(mu_md_up,digits=3),nsmall=3),format(round(mu_md_pi_up,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\n\n\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(exdat$label[5]),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\n#Forest Plot\nif(em==\"SMD\"){\n#SMD\n# \u91cd\u307f\u3092\u683c\u7d0d\u3059\u308b\u30d9\u30af\u30c8\u30eb\u3092\u521d\u671f\u5316\nweights <- numeric(K)\nweight_percentages <- numeric(K)\n\nfor (i in 1:K) {\n  # i\u756a\u76ee\u306e\u7814\u7a76\u306esmd\u306e\u30b5\u30f3\u30d7\u30ea\u30f3\u30b0\u5024\u3092\u53d6\u5f97\n  smd_i_samples <- posterior_samples$delta[, 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_smd_i <- var(smd_i_samples)\n  # \u91cd\u307f\u306f\u5206\u6563\u306e\u9006\u6570\n  weights[i] <- 1 \/ variance_of_smd_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:K,\n  Effective_Weight = weights,\n  Effective_Weight_Percentage = weight_percentages\n)\nprint(results_weights_posterior_var)\nk=K\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\n#setting fs for cex\nfs=1\nif(k>20){fs=round((1-0.02*(k-20)),digits=1)}\n\nm=c(NA,NA,smd,mu_smd,mu_new_smd,NA)\nlw=c(NA,NA,smd_lw,mu_smd_lw,mu_smd_pi_lw,NA)\nup=c(NA,NA,smd_up,mu_smd_up,mu_smd_pi_up,NA)\n\nhete1=\"\"\nhete2=paste(\"I2=\",round(100*mu_smd_I_squared, 2),\"%\",sep=\"\")\n\nhete3=paste(\"tau2=\",format(round(mu_smd_tau_squared,digits=4),nsmall=4),sep=\"\")\nhete4=\"\"\nhete5=paste(\"p=\",p_val_smd, sep=\"\")\nhete6=paste(label,prob_direct_smd, sep=\"\")\n\nrck=rep(NA,K)\nrtk=rck\nfor(i in 1:K)\n{\nrtk[i]=paste(exdat$m1i[i],\" (\",exdat$sd1i[i],\")\",sep=\"\")\nrck[i]=paste(exdat$m2i[i],\" (\",exdat$sd2i[i],\")\",sep=\"\")\n}\n\nau=exdat$author\nsl=c(NA,toString(labe_study),as.vector(au),\"Summary Estimate\",\"Prediction Interval\",NA)\nnc=exdat$n2i\nncl=c(labe_cont,\"Number\",nc,NA,NA,hete1)\nrcl=c(labe_outc,\"Mean (SD)\",rck,NA,NA,hete2)\n\nnt=exdat$n1i\nntl=c(labe_int,\"Number\",nt,NA,NA,hete3)\nrtl=c(labe_outc,\"Mean (SD)\",rtk,NA,NA,hete4)\n\nspac=c(\"    \",NA,rep(NA,k),NA,NA,NA)\nml=c(NA,labe_em,format(round(smd,digits=3),nsmall=3),format(round(mu_smd,digits=3),nsmall=3),NA,hete5)\nll=c(NA,\"95% CI lower\",format(round(smd_lw,digits=3),nsmall=3),format(round(mu_smd_lw,digits=3),nsmall=3),format(round(mu_smd_pi_lw,digits=3),nsmall=3),hete6)\nul=c(NA,\"95% CI upper\",format(round(smd_up,digits=3),nsmall=3),format(round(mu_smd_up,digits=3),nsmall=3),format(round(mu_smd_pi_up,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\n\n\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(exdat$label[5]),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\n##########<\/code><\/pre>\n<\/div><\/details>\n\n\n\n
\"\"<\/figure>\n\n\n\n
metafor\u306efunnel()\u95a2\u6570\u3001regtest()\u3001ranktest()\u95a2\u6570\u3092\u7528\u3044\u3066\u3001Funnel plot\u4f5c\u6210\u3068\u3001Begg\u306e\u691c\u5b9a\u3001Egger\u306e\u691c\u5b9a\u3092\u5b9f\u884c\u3059\u308b\u30b9\u30af\u30ea\u30d7\u30c8\u3092\u4f5c\u6210\u3057\u307e\u3057\u305f\u3002Forest plot\u304c\u51fa\u529b\u3055\u308c\u305f\u5f8c\u3067\u3001\u5b9f\u884c\u3057\u3066\u304f\u3060\u3055\u3044\u3002<\/p>\n\n\n\n
Funnel plot\u4f5c\u6210\u7528\u306e\u30b9\u30af\u30ea\u30d7\u30c8<\/summary>
\n
########################\nlibrary(metafor)\n\n#Plot funnel plot.\nif(em == \"MD\"){\nyi = md\nvi=1\/weights\nsei = sqrt(1\/weights)\nmu = mu_md\ndev.new(width=7,height=7)\nfunnel(yi, sei = sei, level = c(95), refline = mu,xlab = \"MD\", ylab = \"Standard Error\")\n}\nif(em == \"SMD\"){\nyi = smd\nvi=1\/weights\nsei = sqrt(1\/weights)\nmu = mu_smd\ndev.new(width=7,height=7)\nfunnel(yi, sei = sei, level = c(95), refline = mu,xlab = \"SMD\", 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#####################################\n<\/code><\/pre>\n<\/div><\/details>\n\n\n\n
\"\"<\/figure>\n","protected":false},"excerpt":{"rendered":"
\u9023\u7d9a\u5909\u6570\u30a2\u30a6\u30c8\u30ab\u30e0\u3067\u5e73\u5747\u5024\u5deeMean Difference, MD\u307e\u305f\u306f\u6a19\u6e96\u5316\u5e73\u5747\u5024\u5deeStandardized Mean Difference (SMD)\u306e\u7d71\u5408\u5024\u306895\uff05\u78ba\u4fe1\u533a\u9593Credible Interval\u3092\u3001St … “Stan\u3092\u4f7f\u3046\u30d9\u30a4\u30b8\u30a2\u30f3\u30e1\u30bf\u30a2\u30ca\u30ea\u30b7\u30b9\u7528\u306e\u30b3\u30fc\u30c9\uff1a\u9023\u7d9a\u5909\u6570\u30a2\u30a6\u30c8\u30ab\u30e0MD\u3001SMD\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-1893","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\/1893","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=1893"}],"version-history":[{"count":4,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/posts\/1893\/revisions"}],"predecessor-version":[{"id":1911,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=\/wp\/v2\/posts\/1893\/revisions\/1911"}],"wp:attachment":[{"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1893"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1893"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/info.zanet.biz\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1893"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}