Check column names

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"
library(DiagrammeRsvg)
library(rsvg)
## Warning: package 'rsvg' was built under R version 4.4.1
## Linking to librsvg 2.61.0

Load data

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)

Run meta-analysis

res_overall <- rma(yi = EffectSize_num, sei = SE_num, data = dat_clean, method = "REML")
slab_labels <- dat_clean$StudyID

Forest plot (original style restored)

## 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"

Merge effect_size.csv と study_level.csv

# 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

Country から Region を作成

# 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

Region を moderator にしたメタ回帰

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 plot

funnel(res_overall,
       xlab = "Effect size",
       ylab = "Standard Error",
       cex = 0.8)

Egger’s test

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)

Trim and Fill

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
library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)

# 1. Flow 図を SVG として取得
svg <- grViz("
digraph prisma {
  graph [layout = dot, rankdir = TB]
  node [shape = box, style = filled, fillcolor = GhostWhite, fontsize = 12]

  a [label = 'Records identified (n = 2,528)']
  b [label = 'Records after duplicates removed (n = 2,380)']
  c [label = 'Records screened (n = 2,380)']
  d [label = 'Records excluded (n = 2,037)']
  e [label = 'Full-text assessed (n = 343)']
  f [label = 'Full-text excluded (n = 183)']
  g [label = 'Studies meeting criteria (n = 160)']
  h [label = 'Full text available (n = 66)']
  i [label = 'Studies included (n = 16)']
  j [label = 'Effect sizes extracted (k = 55)']

  a -> b
  b -> c
  c -> d
  c -> e
  e -> f
  e -> g
  g -> h
  h -> i
  i -> j
}
") %>% export_svg()

# 2. SVG → PNG に変換
rsvg_png(charToRaw(svg), file = "prisma_flow.png", width = 2000, height = 2000)