1. Data Loading
StudyLevel <- read_csv("Study_Level_0514.csv")
EffectSize <- read_csv("Effect_Size_0514.csv")
StudyLevel <- StudyLevel %>%
rename(Age_Mean = `Age (Mean)`) %>%
mutate(Asia = ifelse(Country == "China", 1, 0))
dat_ls <- EffectSize %>%
left_join(StudyLevel %>% select(StudyID, Age_Mean, Asia, Interval),
by = "StudyID") %>%
mutate(
EffectSize = as.numeric(EffectSize),
SE = as.numeric(SE)
)
dat_ls_clean <- dat_ls %>%
filter(!is.na(EffectSize), !is.na(SE), SE > 0)
df_forest <- dat_ls_clean %>%
mutate(
StudyLabel = paste0("Study_", seq_len(n())),
CI_low = EffectSize - 1.96 * SE,
CI_high = EffectSize + 1.96 * SE
)
2. Main Analysis
res_main <- rma(yi = EffectSize, sei = SE,
data = dat_ls_clean, method = "REML")
summary(res_main)
##
## Random-Effects Model (k = 92; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -66.9020 133.8041 137.8041 142.8258 137.9404
##
## tau^2 (estimated amount of total heterogeneity): 0.2523 (SE = 0.0377)
## tau (square root of estimated tau^2 value): 0.5023
## I^2 (total heterogeneity / total variability): 99.51%
## H^2 (total variability / sampling variability): 205.76
##
## Test for Heterogeneity:
## Q(df = 91) = 10331.6210, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2486 0.0526 4.7253 <.0001 0.1455 0.3517 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
3. Forest Plot（Web + PDF 横向き 1 枚）
3.1 Web（HTML）表示用フォレストプロット
p_web <- ggplot(df_forest, aes(x = EffectSize, y = reorder(StudyLabel, EffectSize))) +
geom_point(size = 2.2, color = "darkblue") +
geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height = 0.15) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Forest Plot (Web View)",
x = "Effect Size (Hedges g)", y = "Study") +
theme_minimal(base_size = 12) +
theme(
axis.text.y = element_text(size = 7, lineheight = 1.1),
plot.margin = margin(15, 15, 15, 15)
)
p_web
3.2 PDF（A4 横向き 1 枚）保存用フォレストプロット
p_pdf <- ggplot(df_forest, aes(x = EffectSize, y = reorder(StudyLabel, EffectSize))) +
geom_point(size = 1.5, color = "darkblue") +
geom_errorbarh(aes(xmin = CI_low, xmax = CI_high), height = 0.08) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
labs(title = "Forest Plot (All Studies, Landscape)",
x = "Effect Size (Hedges g)", y = "Study") +
theme_minimal(base_size = 8) +
theme(
axis.text.y = element_text(size = 5, lineheight = 0.8),
plot.margin = margin(5, 5, 5, 5)
)
ggsave("forest_plot_landscape.pdf", p_pdf,
width = 11.7, height = 8.27)
4. Moderator Analyses（Web + PDF）
4.1 Age
res_age <- rma(yi = EffectSize, sei = SE,
mods = ~ Age_Mean,
data = dat_ls_clean, method = "REML")
summary(res_age)
##
## Mixed-Effects Model (k = 86; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -13.8648 27.7295 57.7295 91.8795 66.3010
##
## tau^2 (estimated amount of residual heterogeneity): 0.0835 (SE = 0.0143)
## tau (square root of estimated tau^2 value): 0.2890
## I^2 (residual heterogeneity / unaccounted variability): 98.60%
## H^2 (unaccounted variability / sampling variability): 71.58
## R^2 (amount of heterogeneity accounted for): 68.31%
##
## Test for Residual Heterogeneity:
## QE(df = 72) = 2526.4365, p-val < .0001
##
## Test of Moderators (coefficients 2:14):
## QM(df = 13) = 190.9487, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.6100 0.1466 4.1597 <.0001 0.3226 0.8974 ***
## Age_Mean12.07 -0.3192 0.1762 -1.8112 0.0701 -0.6647 0.0262 .
## Age_Mean12.56 -0.3742 0.1893 -1.9764 0.0481 -0.7452 -0.0031 *
## Age_Mean12.73 -0.5875 0.2059 -2.8528 0.0043 -0.9911 -0.1839 **
## Age_Mean12.75 -0.4400 0.3265 -1.3475 0.1778 -1.0800 0.2000
## Age_Mean13.15 -0.3814 0.1795 -2.1250 0.0336 -0.7331 -0.0296 *
## Age_Mean14 -0.1841 0.1968 -0.9355 0.3495 -0.5699 0.2016
## Age_Mean15.38 -0.6100 0.1893 -3.2221 0.0013 -0.9811 -0.2389 **
## Age_Mean18.46 -0.6026 0.1793 -3.3617 0.0008 -0.9539 -0.2513 ***
## Age_Mean18.6 -0.3233 0.2240 -1.4434 0.1489 -0.7624 0.1157
## Age_Mean20.5 -0.5893 0.1663 -3.5439 0.0004 -0.9152 -0.2634 ***
## Age_Mean37.9 -0.1360 0.2564 -0.5305 0.5957 -0.6386 0.3665
## Age_Mean61.7 0.8700 0.1796 4.8440 <.0001 0.5180 1.2220 ***
## Age_Mean80+ -0.8259 0.1795 -4.6016 <.0001 -1.1777 -0.4741 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Web
p_age_web <- ggplot(dat_ls_clean, aes(Age_Mean, EffectSize)) +
geom_point(size = 2) +
geom_smooth(method = "lm", color = "red") +
labs(title = "Moderator: Age", x = "Age (Mean)", y = "Effect Size") +
theme_minimal(base_size = 12)
p_age_web
# PDF
p_age_pdf <- ggplot(dat_ls_clean, aes(Age_Mean, EffectSize)) +
geom_point(size = 1.8) +
geom_smooth(method = "lm", color = "red", linewidth = 0.8) +
labs(title = "Moderator: Age", x = "Age (Mean)", y = "Effect Size") +
theme_minimal(base_size = 10)
ggsave("age_moderator.pdf", p_age_pdf,
width = 11.7, height = 8.27)
4.2 Culture
res_culture <- rma(yi = EffectSize, sei = SE,
mods = ~ Asia,
data = dat_ls_clean, method = "REML")
summary(res_culture)
##
## Mixed-Effects Model (k = 86; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -63.4892 126.9784 132.9784 140.2708 133.2784
##
## tau^2 (estimated amount of residual heterogeneity): 0.2630 (SE = 0.0409)
## tau (square root of estimated tau^2 value): 0.5129
## I^2 (residual heterogeneity / unaccounted variability): 99.54%
## H^2 (unaccounted variability / sampling variability): 218.44
## R^2 (amount of heterogeneity accounted for): 0.20%
##
## Test for Residual Heterogeneity:
## QE(df = 84) = 9530.3832, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 1.1631, p-val = 0.2808
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3280 0.0795 4.1246 <.0001 0.1721 0.4838 ***
## Asia -0.1198 0.1111 -1.0785 0.2808 -0.3376 0.0979
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Web
p_culture_web <- ggplot(dat_ls_clean, aes(factor(Asia), EffectSize)) +
geom_boxplot() +
labs(title = "Moderator: Culture (Asia)", x = "Asia (0=Other, 1=China)", y = "Effect Size") +
theme_minimal(base_size = 12)
p_culture_web
# PDF
p_culture_pdf <- ggplot(dat_ls_clean, aes(factor(Asia), EffectSize)) +
geom_boxplot() +
labs(title = "Moderator: Culture (Asia)", x = "Asia (0=Other, 1=China)", y = "Effect Size") +
theme_minimal(base_size = 10)
ggsave("culture_moderator.pdf", p_culture_pdf,
width = 11.7, height = 8.27)
4.3 Interval
res_interval <- rma(yi = EffectSize, sei = SE,
mods = ~ Interval,
data = dat_ls_clean, method = "REML")
summary(res_interval)
##
## Mixed-Effects Model (k = 86; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -17.0330 34.0659 54.0659 77.5040 57.3992
##
## tau^2 (estimated amount of residual heterogeneity): 0.0886 (SE = 0.0147)
## tau (square root of estimated tau^2 value): 0.2977
## I^2 (residual heterogeneity / unaccounted variability): 98.68%
## H^2 (unaccounted variability / sampling variability): 75.55
## R^2 (amount of heterogeneity accounted for): 66.37%
##
## Test for Residual Heterogeneity:
## QE(df = 77) = 2931.2671, p-val < .0001
##
## Test of Moderators (coefficients 2:9):
## QM(df = 8) = 171.1992, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.2908 0.1006 2.8896 0.0039 0.0935 0.4880 **
## Interval10 yr 0.3192 0.1814 1.7596 0.0785 -0.0363 0.6748 .
## Interval12 mo -0.0526 0.1330 -0.3954 0.6925 -0.3133 0.2081
## Interval12 yr 1.1892 0.1467 8.1069 <.0001 0.9017 1.4767 ***
## Interval2 yr -0.5067 0.1465 -3.4577 0.0005 -0.7939 -0.2195 ***
## Interval2.5–3 years -0.2701 0.1290 -2.0938 0.0363 -0.5229 -0.0173 *
## Interval5 mo -0.2908 0.1591 -1.8275 0.0676 -0.6026 0.0211 .
## Interval6 mo -0.1235 0.1220 -1.0125 0.3113 -0.3627 0.1156
## Interval6–9 mo -0.0549 0.1591 -0.3453 0.7299 -0.3668 0.2569
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Web
p_interval_web <- ggplot(dat_ls_clean, aes(Interval, EffectSize)) +
geom_point(size = 2) +
geom_smooth(method = "lm", color = "red") +
labs(title = "Moderator: Interval", x = "Interval (years)", y = "Effect Size") +
theme_minimal(base_size = 12)
p_interval_web
# PDF
p_interval_pdf <- ggplot(dat_ls_clean, aes(Interval, EffectSize)) +
geom_point(size = 1.8) +
geom_smooth(method = "lm", color = "red", linewidth = 0.8) +
labs(title = "Moderator: Interval", x = "Interval (years)", y = "Effect Size") +
theme_minimal(base_size = 10)
ggsave("interval_moderator.pdf", p_interval_pdf,
width = 11.7, height = 8.27)
4.4 Health Category
dat_ls_clean <- dat_ls_clean %>%
mutate(
HealthCat = case_when(
Predictor %in% c("T1 LPAw","T1 LPAh","T1 MVPAw","T1 MVPAh") ~ "PA",
Predictor == "Mobility limitations" ~ "Mobility",
Predictor == "Pain severity" ~ "Pain",
Predictor == "OIDP（poor）" ~ "OralHealth",
Predictor == "SROH（poor）" ~ "SelfRatedHealth",
TRUE ~ "Other"
)
)
res_health_cat2 <- rma(yi = EffectSize, sei = SE,
mods = ~ HealthCat,
data = dat_ls_clean, method = "REML")
summary(res_health_cat2)
##
## Mixed-Effects Model (k = 92; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -52.4120 104.8241 118.8241 136.0045 120.2600
##
## tau^2 (estimated amount of residual heterogeneity): 0.1956 (SE = 0.0302)
## tau (square root of estimated tau^2 value): 0.4423
## I^2 (residual heterogeneity / unaccounted variability): 99.35%
## H^2 (unaccounted variability / sampling variability): 153.13
## R^2 (amount of heterogeneity accounted for): 22.46%
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 6989.1093, p-val < .0001
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 31.0393, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.3100 0.2226 -1.3929 0.1636 -0.7462 0.1262
## HealthCatOralHealth 1.5400 0.3855 3.9951 <.0001 0.7845 2.2955
## HealthCatOther 0.5660 0.2283 2.4790 0.0132 0.1185 1.0135
## HealthCatPA 0.3325 0.3138 1.0597 0.2893 -0.2825 0.9475
## HealthCatPain 0.1875 0.3145 0.5961 0.5511 -0.4289 0.8038
## HealthCatSelfRatedHealth 1.6100 0.3855 4.1767 <.0001 0.8545 2.3655
##
## intrcpt
## HealthCatOralHealth ***
## HealthCatOther *
## HealthCatPA
## HealthCatPain
## HealthCatSelfRatedHealth ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Web
p_health_web <- ggplot(dat_ls_clean, aes(HealthCat, EffectSize)) +
geom_boxplot() +
labs(title = "Moderator: Health Category", x = "Health Category", y = "Effect Size") +
theme_minimal(base_size = 12)
p_health_web
# PDF
p_health_pdf <- ggplot(dat_ls_clean, aes(HealthCat, EffectSize)) +
geom_boxplot() +
labs(title = "Moderator: Health Category", x = "Health Category", y = "Effect Size") +
theme_minimal(base_size = 10)
ggsave("health_category.pdf", p_health_pdf,
width = 11.7, height = 8.27)
library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
library(htmlwidgets)
## 出版バイアスの検討
# Funnel plot
funnel(res_main)
# Egger 検定
regtest(res_main, model = "lm")
##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 3.3282, df = 90, p = 0.0013
## Limit Estimate (as sei -> 0): b = -0.0517 (CI: -0.1797, 0.0763)
# Trim-and-fill
taf <- trimfill(res_main)
taf
##
## Estimated number of missing studies on the left side: 0 (SE = 3.0251)
##
## Random-Effects Model (k = 92; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.2523 (SE = 0.0377)
## tau (square root of estimated tau^2 value): 0.5023
## I^2 (total heterogeneity / total variability): 99.51%
## H^2 (total variability / sampling variability): 205.76
##
## Test for Heterogeneity:
## Q(df = 91) = 10331.6210, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.2486 0.0526 4.7253 <.0001 0.1455 0.3517 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Trim-and-fill 後の Funnel plot
funnel(taf)
# --- ここから追記 ---
# Funnel plot PDF 保存
pdf("funnel_plot.pdf", width = 8, height = 6)
funnel(res_main)
dev.off()
## quartz_off_screen
## 2
# Trim-and-fill 後の Funnel plot PDF 保存
pdf("funnel_plot_trimfill.pdf", width = 8, height = 6)
funnel(taf)
dev.off()
## quartz_off_screen
## 2