dat <- read.csv("effect_size.csv")
names(dat)
## [1] "StudyID" "EffectID" "Outcome.LS.PA.NA."
## [4] "Wave_Comparison" "EffectSize_Type" "EffectSize"
## [7] "SE" "CI_Lower" "CI_Upper"
## [10] "Model_Type" "Adjusted" "Covariates"
## [13] "Age_Effect" "Cohort_Effect" "Period_Effect"
## [16] "Moderator_PA" "Moderator_PA_p" "Notes"
dat <- read.csv("effect_size.csv")
# 数値変換(文字列 → numeric)
dat$EffectSize_num <- as.numeric(dat$EffectSize)
## Warning: NAs introduced by coercion
dat$SE_num <- as.numeric(dat$SE)
## Warning: NAs introduced by coercion
# 変換できなかった行(NA)や SE<=0 を除外
dat_clean <- subset(dat, !is.na(EffectSize_num) & !is.na(SE_num) & SE_num > 0)
res_overall <- rma(yi = EffectSize_num, sei = SE_num, data = dat_clean, method = "REML")
slab_labels <- dat_clean$StudyID
##
Study Level の列名を確認
study <- read.csv("study_level.csv")
names(study)
## [1] "StudyID" "FullCitation"
## [3] "Country" "SampleType"
## [5] "N_Total" "Age_Mean_T1"
## [7] "Age_SD_T1" "Age_Range"
## [9] "Cohort_Info" "Design"
## [11] "Num_Waves" "Interval_Length"
## [13] "LS_Measured" "PA_Measured"
## [15] "NA_Measured" "LS_Scale"
## [17] "PA_Scale" "NA_Scale"
## [19] "PhysicalActivity_Measured" "PA_Type"
## [21] "PA_Coding" "ClinicalSample"
## [23] "Intervention" "Notes"
# Study Level の読み込み
study <- read.csv("study_level.csv")
# effect_size と Study Level を StudyID で結合
dat_merged <- merge(dat_clean, study, by = "StudyID", all.x = TRUE)
# Country の確認
table(dat_merged$Country, useNA = "ifany")
##
## Belgium China Croatia Malaysia UK UK & Germany
## 2 31 1 2 5 3
## USA <NA>
## 7 4
# Asia に含める国のリスト
asia_countries <- c("Japan", "China", "Korea", "South Korea", "Taiwan",
"Singapore", "Hong Kong", "Malaysia", "Thailand")
# Region を作成
dat_merged$Region <- ifelse(dat_merged$Country %in% asia_countries,
"Asia", "NonAsia")
# factor に変換
dat_merged$Region <- factor(dat_merged$Region,
levels = c("NonAsia", "Asia"))
# 確認
table(dat_merged$Region, useNA = "ifany")
##
## NonAsia Asia
## 22 33
res_region <- rma(yi = EffectSize_num,
sei = SE_num,
mods = ~ Region,
data = dat_merged,
method = "REML")
res_region
##
## Mixed-Effects Model (k = 55; tau^2 estimator: REML)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0356 (SE = 0.0075)
## tau (square root of estimated tau^2 value): 0.1886
## I^2 (residual heterogeneity / unaccounted variability): 93.75%
## H^2 (unaccounted variability / sampling variability): 15.99
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 53) = 799.3160, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.2923, p-val = 0.5887
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.0100 0.0423 0.2371 0.8126 -0.0729 0.0930
## RegionAsia 0.0293 0.0543 0.5407 0.5887 -0.0770 0.1357
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
funnel(res_overall,
xlab = "Effect size",
ylab = "Standard Error",
cex = 0.8)
egger_test <- regtest(res_overall, model = "rma")
egger_test
##
## Regression Test for Funnel Plot Asymmetry
##
## Model: mixed-effects meta-regression model
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: z = -2.4797, p = 0.0131
## Limit Estimate (as sei -> 0): b = 0.1773 (CI: 0.0486, 0.3059)
tf <- trimfill(res_overall)
tf
##
## Estimated number of missing studies on the right side: 0 (SE = 4.3553)
##
## Random-Effects Model (k = 55; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.0349 (SE = 0.0073)
## tau (square root of estimated tau^2 value): 0.1869
## I^2 (total heterogeneity / total variability): 93.65%
## H^2 (total variability / sampling variability): 15.74
##
## Test for Heterogeneity:
## Q(df = 54) = 799.5702, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0279 0.0263 1.0636 0.2875 -0.0235 0.0794
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1