#3 practice
Sangyoung
2023-08-07
1 package loading
2 data loading
dat <- read_excel("Practice Dataset Volume 07082023 (1).xlsx") %>%
mutate(region_wb = countrycode(country, origin = 'country.name', destination = c('region'))) %>% # regional data from WB
inner_join(wbstats::wb_cachelist$countries %>%
select(country, iso3c, income_level),by = "country")3 Exploratory data Analysis
3.1 Word cloud
interventions at a glace (focusing on frequently used words in intervention description)
library(tidytext)
library(SnowballC)
dat %>%
tidytext::unnest_tokens(word, inter_desc) %>%
anti_join(stop_words) %>%
mutate(word = wordStem(word)) %>%
filter(str_detect(word, "[::digits::]")) %>%
count(word, sort = T) %>%
filter(n > 10) %>%
wordcloud2::wordcloud2()## Joining with `by = join_by(word)`- Based on the word cloud analysis to gauge the interventions, “Road” appears to be a predominant term, suggesting that the intervention might pertain to road construction. Additionally, words like “infrastructure” and “bridge” also emerged from this analysis.
3.2 Summary of studies
library(plotly)
a <- dat %>%
mutate(authors_rsvd = str_remove_all(authors, "\\s+\\d+")) %>%
mutate(authors_rsvd = str_replace_all(authors_rsvd, "and", "&")) %>%
mutate(authors_year = str_extract(authors, "\\d+")) %>%
mutate(authors_rsvd = glue::glue("{authors_rsvd} ({authors_year})")) %>%
# select(yi, nfci_lower, nfci_upper, n_T_revised, authors_ctry_rob_inter) %>%
mutate(authors_rsvd = fct_reorder(authors_rsvd , yi)) %>%
ggplot(aes(yi, authors_rsvd,
text = paste0(
"Intervention subcode =",inter_subcode,"\n",
"ROB rating =", overall_rob,"\n",
"Region=", region,"\n",
"Country =", country,"\n",
"Outcome Description =", outcome_desc,"\n",
"Sample size =", n_T_revised))) +
geom_errorbarh(aes(xmin = nfci_lower,
xmax = nfci_upper),
height = 0.3,
size = 2) +
geom_point(size = 3, show.legend = F) +
geom_vline(xintercept = 0, lty = 2, size = 2, , show.legend = F, color = "red", alpha = .7) +
labs(y = "") +
labs(x = "Weighted average effect and confidence interval") +
theme(axis.text.y = element_text(size = 20, hjust = 0))
b <- ggplotly(a, tooltip = c("text"))
layout(b, yaxis = list(categoryorder = "array"))- This is a concise version of a forest plot with additional details. This interactive plot allows for the identification of key elements such as “intervention subgroup,” “ROB rating,” “Region,” “Country,” “Outcome description,” and “Sample size.”
- Confidence interval in this plot represented “nfci_upper” and “nfci_lower”, which calculated yi +_ 1.96* sqrt(vi)
3.3 data check
3.3.1 Intervention sub group
dat %>%
count(inter_subcode) %>%
mutate(pct = n / 16) %>%
mutate(inter_subcode = fct_reorder(inter_subcode, pct)) %>%
ggplot(aes(pct, glue::glue("{inter_subcode} \n(N = {n})"), fill = inter_subcode)) +
geom_col(show.legend = F) +
geom_label(aes(label = scales::percent(pct)), show.legend = F) +
labs(x = "Percentage",
y = "") +
theme(axis.text.y = element_text(hjust = 0.5))
dat %>%
filter(inter_subcode == "Arranged or curated offers from buyers") %>%
pull(inter_desc)## [1] "p.1: The SDCP, implemented by the Government of Kenya (GoK) from 2005 to 2015, was\r\ndesigned to reach dairy groups engaged in milk production in nine milk-producing\r\ncounties. The overall goal of SDCP was to increase the income of poor rural households\r\nthat depend substantially on production and trade of dairy products for their livelihoods.\r\nThe SDCP was implemented through various interrelated components.\r\nFirst, SDCP provided training on organisational, managerial, and enterprise skills (e.g.\r\nbookkeeping, accounting and financial planning) to farmers. Second, the programme\r\ntargeted household production by aiming to enhance dairy farming productivity and\r\nreduce production costs through grants, trainings, field days and demonstrations. Last,\r\nSDCP aimed to strengthen market linkages for small-scale milk producers by improving\r\nroad infrastructure and conducting additional trainings on milk handling practices and\r\nvalue-addition opportunities."- Not all interventions pertain to “road construction.” Another intervention involves “Facilitated or curated offers from buyers,” which appears to influence milk production.
3.3.2 Intervention sub group
dat %>%
select(multi_inter, inter_subcode, id)## # A tibble: 16 × 3
## multi_inter inter_subcode id
## <dbl> <chr> <chr>
## 1 0 Domestic transport infrastructure 7362…
## 2 0 Domestic transport infrastructure 7362…
## 3 0 Domestic transport infrastructure 7362…
## 4 1 Arranged or curated offers from buyers 7362…
## 5 0 Domestic transport infrastructure 7362…
## 6 0 Domestic transport infrastructure 7362…
## 7 1 Domestic transport infrastructure 7379…
## 8 0 Domestic transport infrastructure 7380…
## 9 0 Domestic transport infrastructure 7380…
## 10 1 Domestic transport infrastructure 7466…
## 11 1 Domestic transport infrastructure ; Alternative physical … 7471…
## 12 0 Domestic transport infrastructure 7471…
## 13 1 Domestic transport infrastructure 7471…
## 14 0 Domestic transport infrastructure 7471…
## 15 1 Domestic transport infrastructure 7541…
## 16 0 Domestic transport infrastructure 7727…- There appears to be an inconsistency in the “multi_inter” column. Entries such as “73797268_4S8S12S16” and “74666345_19” are coded as “1”, indicating multiple interventions, even though only one intervention subgroup is evident.
3.3.3 Income level using WB category
dat %>%
count(income_level) %>%
mutate(pct = n / 16) %>%
mutate(income_level = fct_reorder(income_level, pct)) %>%
ggplot(aes(pct, income_level, fill = income_level)) +
geom_col(show.legend = F) +
geom_label(aes(label = scales::percent(pct)), show.legend = F) +
labs(x = "Percentage",
y = "") +
theme(axis.text.y = element_text(hjust = 0.5))
- One country (Trinidad and Tobago) comes from high income country
“The twin island state of Trinidad and Tobago, is a small but high-income country. With a population of 1.3 million people and rich in natural resources, it has one of the highest GNI per capita in Latin America and the Caribbean (US$17,640 in 2015, Atlas method)” (WB website) https://www.worldbank.org/en/country/trinidadandtobago/overview#:~:text=The%20twin%20island%20state%20of,in%202015%2C%20Atlas%20method).
3.3.4 Countries and regions
dat %>%
count(country, region, region_wb)## # A tibble: 15 × 4
## country region region_wb n
## <chr> <chr> <chr> <int>
## 1 Argentina Latin America Latin Americ… 1
## 2 Armenia East Europe and Central and West Asia Europe & Cen… 1
## 3 Bangladesh South-East Asia South Asia 2
## 4 China East Asia East Asia & … 1
## 5 El Salvador Latin America Latin Americ… 1
## 6 India South Asia South Asia 1
## 7 Kenya Sub-Saharan Africa Sub-Saharan … 1
## 8 Liberia Sub-Saharan Africa Sub-Saharan … 1
## 9 Madagascar Sub-Saharan Africa Sub-Saharan … 1
## 10 Nicaragua Latin America Latin Americ… 1
## 11 Nigeria Sub-Saharan Africa Sub-Saharan … 1
## 12 Rwanda Sub-Saharan Africa Sub-Saharan … 1
## 13 Tanzania Sub-Saharan Africa Sub-Saharan … 1
## 14 Trinidad and Tobago Latin America Latin Americ… 1
## 15 Zambia Sub-Saharan Africa Sub-Saharan … 1- Continent categories seem differ between WB category and the coded region. It would be great to specify the regional cateogires. For instance, Bangladesh was coded “South-East Asia” in the previous coding, whereas the WB put it to “South Asia”.
3.4 correlation plots to show the correlation between yi (effect size) and other variables. It might help figure out some moderators
3.4.1 Continuous variable
library(corrr)## Warning: 패키지 'corrr'는 R 버전 4.2.3에서 작성되었습니다dat %>%
select_if(is.numeric) %>%
select(-multi_inter) %>%
correlate() %>%
focus(yi) %>%
mutate(direction = ifelse(yi > 0, "Positive", "Negative")) %>%
mutate(term = fct_reorder(factor(term), yi)) %>%
ggplot(aes(yi, term, fill =direction)) +
geom_col() +
labs(x = "Correlation with effect size (yi)",
y= "") +
theme(axis.text.y = element_text(hjust = .5))## Warning in stats::cor(x = x, y = y, use = use, method = method): 표준편차가
## 0입니다## Correlation computed with
## • Method: 'pearson'
## • Missing treated using: 'pairwise.complete.obs'## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `term = fct_reorder(factor(term), yi)`.
## Caused by warning:
## ! `fct_reorder()` removing 1 missing value.
## ℹ Use `.na_rm = TRUE` to silence this message.
## ℹ Use `.na_rm = FALSE` to preserve NAs.## Warning: Removed 1 rows containing missing values (`position_stack()`).
- There don’t appear to be noticeable variables correlated with the effect size (yi), aside from certain statistical methods used to determine the effect size.
4 Modeling
4.1 Road intervention (k = 16)
- The estimates show positive effects but are not statistically significant with a p-value of 0.4344.
4.2 forest plot
#forest plot ordered by effect size
resoutcomesubcode_meta <- metagen(TE = yi,
seTE = sqrt(vi),
data=dat,
sm = "SMD",
fixed = F,
random = T,
method.tau = "DL",
studlab = authors_ctry_rob_inter)
forest.meta(resoutcomesubcode_meta,
sortvar = TE,
leftcols = c("studlab"),
leftlabs = c("Author (Year), Country, Intervention abbreviation, Crop "),
rightcols = c("TE", "ci"),
weight.study = "common",
digits.TE = 2,
print.Q = F,
print.I2 = F,
print.pval.Q = F,
print.tau2 = F,
random =F,
text.addline2 = format("\n\n* indicates some concerns or low risk of bias"))
- Although there are five studies excluding “high risk of bias”, but this plot only reports four studies indicating “some concern” and “low risk of bias”
4.3 Outlier assessment
inf <- influence(resoutcomesubcode)
plot(inf)
- There is no indication of outliers
4.4 Publication bias
funnel(resoutcomesubcode)
funnel(trimfill(resoutcomesubcode))
- trimfill results indicate two missing outcomes. Two unfilled dots on the left represented the possible missing estimates.
regtest(resoutcomesubcode, model = "lm")##
## Regression Test for Funnel Plot Asymmetry
##
## Model: weighted regression with multiplicative dispersion
## Predictor: standard error
##
## Test for Funnel Plot Asymmetry: t = 0.6772, df = 14, p = 0.5093
## Limit Estimate (as sei -> 0): b = -0.0226 (CI: -0.1589, 0.1137)- With a p-value grater than 0.05, we fail to reject the hypothesis that the funnel plot is asymmetric. This suggests the presence of publication bias.
- ” To evaluate the funnel asymmetry, we inspect the size of B0, and if it differs significantly from zero. When this is the case, Egger’s test indicates funnel plot asymmetry.
https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/pub-bias.html?#eggers-test
4.5 moderator analysis
(previously, the command did not specify the method. But it would be great to use the same estimator as the main regression, add the prompt “method =”DL”)
4.5.1 Region using World bank category (not significant)
rma(yi, vi,
mods = ~ region_wb,
method = "DL", # the same method as the main meta regression
data = dat %>%
mutate(region_wb = fct_relevel(region_wb, "Sub-Saharan Africa")) # assign the reference group in the regression
) %>%
tidy()## # A tibble: 5 × 6
## term type estimate std.error statis…¹ p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 intercept summary 0.0394 0.0473 0.834 0.404
## 2 region_wbEast Asia & Pacific summary -0.199 0.123 -1.61 0.107
## 3 region_wbEurope & Central Asia summary -0.0325 0.118 -0.276 0.782
## 4 region_wbLatin America & Caribbean summary 0.0519 0.0820 0.633 0.527
## 5 region_wbSouth Asia summary -0.0546 0.0891 -0.613 0.540
## # … with abbreviated variable name ¹statisticdat %>%
mutate(region_wb = fct_relevel(region_wb, "Sub-Saharan Africa")) %>%
count(region_wb)## # A tibble: 5 × 2
## region_wb n
## <fct> <int>
## 1 Sub-Saharan Africa 7
## 2 East Asia & Pacific 1
## 3 Europe & Central Asia 1
## 4 Latin America & Caribbean 4
## 5 South Asia 3- reference group (Sub-Saharan Africa (N = 7, the largest)). there is no statistical significance.
4.5.2 income level category by WB (not significant)
rma(yi, vi,
mods = ~ income_level,
method = "DL", # the same method as the main meta regression
data = dat %>%
mutate(income_level = fct_relevel(income_level, c("Low income", "Lower middle income", "Upper middle income", "High income"))))##
## Mixed-Effects Model (k = 16; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0101 (SE = 0.0060)
## tau (square root of estimated tau^2 value): 0.1003
## I^2 (residual heterogeneity / unaccounted variability): 75.65%
## H^2 (unaccounted variability / sampling variability): 4.11
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 12) = 49.2851, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 2.7445, p-val = 0.4327
##
## Model Results:
##
## estimate se zval pval ci.lb
## intrcpt 0.0440 0.0660 0.6668 0.5049 -0.0854
## income_levelLower middle income -0.0318 0.0779 -0.4085 0.6829 -0.1845
## income_levelUpper middle income -0.0722 0.0984 -0.7337 0.4631 -0.2650
## income_levelHigh income 0.1635 0.1439 1.1364 0.2558 -0.1185
## ci.ub
## intrcpt 0.1734
## income_levelLower middle income 0.1208
## income_levelUpper middle income 0.1207
## income_levelHigh income 0.4455
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dat %>%
count(income_level)## # A tibble: 4 × 2
## income_level n
## <chr> <int>
## 1 High income 1
## 2 Low income 3
## 3 Lower middle income 9
## 4 Upper middle income 3- reference group (Low income (N = 3)). there is no statistical significance.
- income level might be associated with the starting point of intervention. Since the intervention is road construction, we might assume the country wealth can be associated with the initial infrastructure.
4.5.3 Multiple interventions (Not significant)
rma(yi, vi,
mods = ~ multi_inter2,
method = "DL", # the same method as the main meta regression
data = dat %>%
mutate(multi_inter2 = ifelse(str_detect(inter_subcode,";"), "No", "Yes"))
)##
## Mixed-Effects Model (k = 16; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0065 (SE = 0.0041)
## tau (square root of estimated tau^2 value): 0.0805
## I^2 (residual heterogeneity / unaccounted variability): 68.71%
## H^2 (unaccounted variability / sampling variability): 3.20
## R^2 (amount of heterogeneity accounted for): 24.81%
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 44.7488, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 4.1061, p-val = 0.0427
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt -0.1597 0.0933 -1.7125 0.0868 -0.3425 0.0231 .
## multi_inter2Yes 0.1969 0.0971 2.0264 0.0427 0.0064 0.3873 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dat %>%
mutate(multi_inter2 = ifelse(str_detect(inter_subcode,";"), "No", "Yes")) %>%
count(multi_inter2)## # A tibble: 2 × 2
## multi_inter2 n
## <chr> <int>
## 1 No 1
## 2 Yes 15- reference group (Single intervention (N = 15))
- Multi-intervention can offer certain advantages, particularly in combining hardware and software. However, its statistical significance is at the 0.05 level, although this observation may be influenced by a single study.
4.5.4 Road and Bridge
4.5.4.1 Road (Not significant)
rma(yi, vi,
mods = ~ road,
method = "DL", # the same method as the main meta regression
data = dat %>%
mutate(road = str_detect(road_bridge, "Road"))
)##
## Mixed-Effects Model (k = 16; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0100 (SE = 0.0057)
## tau (square root of estimated tau^2 value): 0.0998
## I^2 (residual heterogeneity / unaccounted variability): 76.61%
## H^2 (unaccounted variability / sampling variability): 4.28
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 59.8664, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.0018, p-val = 0.9666
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.0261 0.0773 0.3378 0.7355 -0.1254 0.1776
## roadTRUE -0.0035 0.0841 -0.0419 0.9666 -0.1684 0.1614
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dat %>%
mutate(road = str_detect(road_bridge, "Road")) %>% count(road)## # A tibble: 2 × 2
## road n
## <lgl> <int>
## 1 FALSE 2
## 2 TRUE 14- reference group (No road intervention (N = 2, Bridge))
4.5.4.2 Bridge and road (Not significant)
rma(yi, vi,
mods = ~ road_bridge,
method = "DL", # the same method as the main meta regression
data = dat %>%
mutate(road_bridge = str_detect(road_bridge, ";"))
)##
## Mixed-Effects Model (k = 16; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0087 (SE = 0.0051)
## tau (square root of estimated tau^2 value): 0.0935
## I^2 (residual heterogeneity / unaccounted variability): 75.15%
## H^2 (unaccounted variability / sampling variability): 4.02
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 56.3458, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.9255, p-val = 0.3360
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.0056 0.0341 0.1630 0.8705 -0.0613 0.0724
## road_bridgeTRUE 0.0629 0.0654 0.9621 0.3360 -0.0652 0.1910
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- reference group (No road intervention (N = 12, Bridge))
- Bridge and road construction can be associated with the synergy, but it is not significant
4.5.5 COVID 19 effect?
- whether intervention period include 2019 or not
rma(yi, vi,
mods = ~ year_end,
method = "DL", # the same method as the main meta regression
data = dat %>%
mutate(year_end = year_inter + length_followup/12 + expfollow_up) %>%
mutate(year_end = year_end >= 2019))## Warning: 1 study with NAs omitted from model fitting.##
## Mixed-Effects Model (k = 15; tau^2 estimator: DL)
##
## tau^2 (estimated amount of residual heterogeneity): 0.0095 (SE = 0.0055)
## tau (square root of estimated tau^2 value): 0.0972
## I^2 (residual heterogeneity / unaccounted variability): 76.77%
## H^2 (unaccounted variability / sampling variability): 4.31
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 13) = 55.9668, p-val < .0001
##
## Test of Moderators (coefficient 2):
## QM(df = 1) = 0.2882, p-val = 0.5914
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.0131 0.0329 0.3997 0.6894 -0.0513 0.0776
## year_endTRUE 0.0434 0.0809 0.5369 0.5914 -0.1151 0.2019
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dat %>%
mutate(year_end = year_inter + length_followup/12 + expfollow_up) %>%
mutate(year_end = year_end >= 2019) %>%
count(year_end)## # A tibble: 3 × 2
## year_end n
## <lgl> <int>
## 1 FALSE 13
## 2 TRUE 2
## 3 NA 1- “73628861_22” need to revist the data entry, because the end year is 2027
- Reference groups (intervention ended before 2019, N = 13) _ Despite the potential negative effects of the COVID-19 pandemic on outcomes, it appears that in this case, the outcomes have not been significantly impacted.
5 Risk of bias
datlowrob <- dat %>%
filter(!str_detect(overall_rob, "High"))resoutcomesubcodenohi <- rma(yi, vi, method="DL", data=datlowrob)
resoutcomesubcodenohi##
## Random-Effects Model (k = 5; tau^2 estimator: DL)
##
## tau^2 (estimated amount of total heterogeneity): 0.0226 (SE = 0.0188)
## tau (square root of estimated tau^2 value): 0.1502
## I^2 (total heterogeneity / total variability): 87.22%
## H^2 (total variability / sampling variability): 7.83
##
## Test for Heterogeneity:
## Q(df = 4) = 31.3010, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.0485 0.0723 0.6702 0.5027 -0.0932 0.1902
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1resoutcomesubcodehi <- rma(yi, vi, method="DL", data= dat %>% filter(str_detect(overall_rob, "High")))
resoutcomesubcodenohi_meta <- metagen(TE = yi,
seTE = sqrt(vi),
data=datlowrob,
sm = "SMD",
fixed = F,
random = T,
method.tau = "DL",
studlab = authors_ctry_rob_inter)
forest.meta(resoutcomesubcodenohi_meta,
sortvar = TE,
leftcols = c("studlab"),
leftlabs = c("Author (Year), Country, Intervention abbreviation, Crop "),
# rightcols = c("TE", "ci", "w.random"),
weight.study = "common",
digits.TE = 2,
print.Q = T,
print.I2 = T,
print.pval.Q = T,
text.addline2 = format("\n\n* indicates some concerns or low risk of bias"))
Number_study <-
dat %>%
mutate(type = ifelse(overall_rob != "High risk of bias", "estimates \nexcluding high risk of bias", "others")) %>%
count(type) %>%
janitor::adorn_totals(name = "all estimates") %>%
as_tibble() %>%
filter(type != "others")
comparision <-
tidy(resoutcomesubcode, conf.int = T) %>%
mutate(type = "all estimates") %>%
# bind_rows(tidy(resoutcomesubcodehi, conf.int = T) %>%
# mutate(type = "estimates only for high risk of bias"))
bind_rows(tidy(resoutcomesubcodenohi, conf.int = T) %>%
mutate(type = "estimates \nexcluding high risk of bias")) %>%
inner_join(Number_study, by = "type")
comparision %>%
mutate(type = str_to_sentence(type)) %>%
mutate(type = glue::glue("{type} (N = {n})")) %>%
mutate(type = fct_reorder(type, estimate)) %>%
ggplot(aes(x = estimate , y = type, color = type)) +
geom_errorbarh(aes(xmin = conf.low,
xmax = conf.high),
height = 0.3,
size = 2, show.legend = F) +
geom_point(size = 5, show.legend = F) +
geom_vline(xintercept = 0, lty = 2, linewidth = 2, , show.legend = F) +
labs(y = "",
x = "Weighted average effect and confidence interval") +
theme_light() +
theme(axis.text.y = element_text(size = 20,
hjust = .5)) 
- Excluding the high-risk study leads to slightly higher estimates. The confidence interval appears wider compared to the others due to the reduction in sample size. To summarize, both estimates encompass 0 within their confidence intervals, indicating that there is no significant change in the statistical significance.
6 Some points to discuss
6.1 Different random effect methods (DL vs REML)
- Estimate seems different, so it would be great to use the same model specification (method = “DL”)
6.2 Reference group
rma(yi, vi,
mods = ~ nonstaggered,
data=dat) %>% tidy()## # A tibble: 3 × 6
## term type estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 intercept summary 0.503 0.408 1.23 0.218
## 2 nonstaggeredNo summary -0.484 0.410 -1.18 0.238
## 3 nonstaggeredYes summary -0.475 0.414 -1.15 0.251dat %>% count(nonstaggered)## # A tibble: 3 × 2
## nonstaggered n
## <chr> <int>
## 1 NA 1
## 2 No 12
## 3 Yes 3- R took “NA” as a string outcome, so the result indicated the three outcomes (NA, nonstaggeredNo, nonstaggeredYes)
rma(yi, vi,
mods = ~ nonstaggered,
data=dat %>% filter(nonstaggered != "NA")) %>% tidy()## # A tibble: 2 × 6
## term type estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 intercept summary 0.0188 0.0371 0.506 0.613
## 2 nonstaggeredYes summary 0.00833 0.0775 0.108 0.914If we filtered out “NA” observations, the results can be more intepretable. non-staggeredYes vs non-staggeredNo
So, it would be better if the reference groups need to be reported in the moderator analysis