fp_outcome_meta_analysis
James McLauchlan
2025-07-04
citation
## Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software, 36(3), 1-48. doi:10.18637/jss.v036.i03
## viz funnel url: https://cloud.r-project.org/web/packages/metaviz/vignettes/metaviz.htmlViolent reoffending
#load excel file
forensic_patient_violent<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 2)
#View(forensic_patient_violent)
## glimpse(forensic_patient_violent)
#str(forensic_patient_violent)
forensic_patient_violent <- forensic_patient_violent %>%
rename(after_2000 = after_2000)
forensic_patient_violent <- forensic_patient_violent %>%
rename(secure_level = `security level`)random effect model for violent reoffending
forensic_patient_meta_result_violent <- metarate(event = violent_reoffend,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_violent,
sm = "IRLN", # Log incidence rate
method = "Inverse", # DerSimonian-Laird (default) for random effects
common = FALSE, # Use random-effects model
irscale = 100000, # per 100,000 person-years
random = TRUE)
summary(forensic_patient_meta_result_violent)## events 95%-CI %W(random)
## Coid _2007 3618.5174 [ 3194.2582; 4099.1266] 5.7
## Bengtson_2019 2380.9524 [ 2046.8173; 2769.6337] 5.7
## Baxter_1999 16687.0167 [12286.9291; 22662.8251] 5.6
## Dean_2021* 144.6655 [ 54.2955; 385.4478] 4.8
## Blattner_2009 3139.7174 [ 1410.5516; 6988.6314] 5.0
## Fazel_2016 2574.9424 [ 2478.0819; 2675.5888] 5.7
## Thomson_2023 413.9073 [ 235.0620; 728.8259] 5.3
## Simpson_2006 6944.4444 [ 3119.8659; 15457.4939] 5.0
## Dowset_2005 6808.5106 [ 3404.9201; 13614.3626] 5.2
## Hogan_2019 2082.0940 [ 1233.1251; 3515.5519] 5.4
## deVogel_2019 1521.0778 [ 900.8620; 2568.2932] 5.4
## Miraglia_2011 740.1925 [ 542.9475; 1009.0936] 5.6
## Sivak_2023 3993.0556 [ 3153.7897; 5055.6614] 5.6
## Dolan_2004 1863.3540 [ 837.1317; 4147.6008] 5.0
## Simpson_2018 16666.6667 [ 8967.5785; 30975.7843] 5.3
## Ojansuu_2023 1165.3314 [ 878.1912; 1546.3571] 5.6
## Edwards_2002 303.0303 [ 42.6859; 2151.2338] 3.2
## Davison_1999 5000.0000 [ 3063.1596; 8161.5076] 5.4
## Friendship_1999 2763.4839 [ 1920.4048; 3976.6841] 5.5
##
## Number of studies: k = 19
## Number of events: e = 3352
##
## events 95%-CI
## Random effects model 2375.9045 [1399.5739; 4033.3148]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 1.2817 [0.7046; 3.1926]; tau = 1.1321 [0.8394; 1.7868]
## I^2 = 95.7% [94.3%; 96.7%]; H = 4.80 [4.20; 5.49]
##
## Test of heterogeneity:
## Q d.f. p-value
## 414.86 18 < 0.0001
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Log transformation
## - Normal approximation confidence interval for individual studies
## - Events per 100000 person-years## pooled rate = 2375 per 100,000 person-years
## metareg(forensic_patient_meta_result_violent, ~ total_year_complete)Plots
forest(forensic_patient_meta_result_violent, backtransf = TRUE,
xlab = "Incidence rate per 100,000 person-years",
refline = exp(forensic_patient_meta_result_violent$b),
addfit = TRUE,
xlim = c(-200, 35000))
##These values suggest large variation in the true incidence rates across studies, even after accounting for within-study error
funnel(forensic_patient_meta_result_violent,
xlim = c(-10, 10))
funnel(forensic_patient_meta_result_violent,
xlim = c(-8, 8),
xlab = "Log Incidence Rate",
contour = c(0.90, 0.95, 0.99),
col.contour = c("gray80", "gray60", "gray40"),
backtransf = FALSE)
legend("topright", # position
legend = c("p < 0.01", "p < 0.05", "p < 0.10"),
fill = c("gray40", "gray60", "gray80"),
border = TRUE,
bty = "o",
title = "Significance regions",
cex = 0.9)
#names(forensic_patient_meta_result_violent)
#str(forensic_patient_meta_result_violent)transformation for meta regression - violent reoffending
## Moderators: mean age (<37.5 vs.>37.5), ssd percentage (< 71 vs. >71), UK studies (UK vs. non-UK) security level (0-2), publication year, country, , , follow-up duration
print(sort(forensic_patient_violent$m_age))## [1] 31.6 32.7 33.6 33.8 36.4 36.6 37.4 37.5 37.7 37.7 38.6 41.5 42.3 44.0 44.5
## [16] 47.0forensic_patient_violent <- forensic_patient_violent %>%
mutate(age_moderator = if_else(m_age < 37.5, 0, 1))
forensic_patient_violent <- forensic_patient_violent %>%
mutate(age_moderator = case_when(
!is.na(m_age) & m_age < 37.5 ~ 0,
!is.na(m_age) & m_age >= 37.5 ~ 1,
TRUE ~ NA_real_
))
#table(forensic_patient_violent$age_moderator)
forensic_patient_violent <- forensic_patient_violent %>%
mutate(ssd_percent = if_else(!is.na(ssd) & !is.na(sample), #checks for fields with no numbers
(ssd / sample) * 100,
NA_real_))
#print(sort(forensic_patient_violent$ssd_percent))
forensic_patient_violent <- forensic_patient_violent %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 71, 1, 0)))
forensic_patient_violent <- forensic_patient_violent %>%
mutate(age_moderator = if_else(m_age < 37.5, 0, 1))
forensic_patient_violent <- forensic_patient_violent %>%
mutate(ssd_percent = if_else(!is.na(ssd) & !is.na(sample), #checks for fields with no numbers
(ssd / sample) * 100,
NA_real_))
#print(sort(forensic_patient_violent$ssd_percent))
forensic_patient_violent <- forensic_patient_violent %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 71, 1, 0)))
#table(forensic_patient_violent$ssd_high, useNA = "ifany")
# forensic_patient_violent <- forensic_patient_violent %>%
#mutate(men_percent = if_else(!is.na(men) & !is.na(sample),
#(men / sample) * 100,
# NA_real_))
##print(sort(forensic_patient_violent$men_percent)) > abandoned as the split would not be meaningful. Most studies all but 3 studies had less than 80% men
forensic_patient_violent <- forensic_patient_violent %>%
mutate(uk_studies = if_else(tolower(country) == "uk", 1, 0))
#table(forensic_patient_violent$uk_studies, useNA = "ifany")Meta-regression - violent reoffending
forensic_patient_violent <- forensic_patient_violent %>%
filter(!is.na(violent_reoffend), !is.na(total_year_complete),
violent_reoffend > 0, total_year_complete > 0) %>%
mutate(
rate = violent_reoffend / total_year_complete * 100000, # per 100,000 person-years
yi = log(rate),
vi = 1 / violent_reoffend # Poisson approximation for log incidence rate variance
)
# Step 3: Run meta-regression
res <- rma(yi = yi, vi = vi,
mods = ~ secure_level + after_2000 + age_moderator + ssd_high + uk_studies + m_follow,
data = forensic_patient_violent,
method = "REML")## Warning: 9 studies with NAs omitted from model fitting.## Warning: Redundant predictors dropped from the model.summary(res)##
## Mixed-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -3.7162 7.4323 21.4323 17.1364 133.4323
##
## tau^2 (estimated amount of residual heterogeneity): 0.3259 (SE = 0.2638)
## tau (square root of estimated tau^2 value): 0.5709
## I^2 (residual heterogeneity / unaccounted variability): 91.50%
## H^2 (unaccounted variability / sampling variability): 11.77
## R^2 (amount of heterogeneity accounted for): 67.76%
##
## Test for Residual Heterogeneity:
## QE(df = 4) = 32.6723, p-val < .0001
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 20.7545, p-val = 0.0009
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 10.6141 0.9188 11.5520 <.0001 8.8133 12.4150 ***
## secure_level 0.0176 0.2043 0.0864 0.9312 -0.3828 0.4181
## age_moderator -0.7803 0.8672 -0.8998 0.3682 -2.4800 0.9194
## ssd_high -0.6776 0.6421 -1.0553 0.2913 -1.9361 0.5809
## uk_studies -0.8859 0.5116 -1.7315 0.0834 -1.8886 0.1169 .
## m_follow -0.1753 0.0529 -3.3155 0.0009 -0.2789 -0.0717 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##colnames(model.matrix(res))reoffending
#load excel file
forensic_patient_reoffend<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 3)
#View(forensic_patient_reoffend)
#glimpse(forensic_patient_reoffend)
forensic_patient_reoffend <- forensic_patient_reoffend %>%
rename(secure_level = `security level(0=unspecified or mixed/1=low/2=medium/3=high)`)
forensic_patient_reoffend <- forensic_patient_reoffend %>%
rename(after_2000 = `after_2000 (0=no/1=yes)`)random effect model for reoffending
forensic_patient_meta_result_reoffend <- metarate(event = reoffending,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_reoffend,
sm = "IRLN", # Log incidence rate
method = "Inverse", # DerSimonian-Laird (default) for random effects
common = FALSE, # Use random-effects model,
irscale = 100000, # per 100,000 person-years
random = TRUE)
summary(forensic_patient_meta_result_reoffend)## events 95%-CI %W(random)
## Argo_2023 210.0840 [ 52.5415; 840.0080] 3.9
## Fioritti_2011 2836.8794 [1064.7315; 7558.6050] 4.9
## Edwards_2002 2121.2121 [1011.2538; 4449.4675] 5.5
## Nilsson_2011 4227.0531 [2015.1796; 8866.6926] 5.5
## Dowset_2005 6808.5106 [3404.9201; 13614.3626] 5.6
## Simpson_2006 10416.6667 [5419.9455; 20019.9328] 5.7
## Halstead_2001 9243.6975 [5119.1612; 16691.3954] 5.8
## Probst_2020 5476.4513 [3301.5656; 9084.0292] 6.0
## Dean_2021* 614.8282 [ 382.2144; 989.0095] 6.1
## Russo_1994 4395.6044 [2835.8558; 6813.2300] 6.2
## Philipse_2006 3176.9306 [2163.0839; 4665.9715] 6.2
## Tellefsen_1992 10666.6667 [7543.2088; 15083.4719] 6.3
## Marshall_2014 4494.3820 [3386.9565; 5963.8999] 6.4
## Seifert_2005 4235.0042 [3251.4558; 5516.0709] 6.4
## Silver_1989 5174.3532 [4086.8006; 6551.3182] 6.5
## Miraglia_2011 1480.3849 [1189.0700; 1843.0702] 6.5
## Sivak_2023 7638.8889 [6440.8414; 9059.7826] 6.5
##
## Number of studies: k = 17
## Number of events: e = 542
##
## events 95%-CI
## Random effects model 3744.7390 [2443.3231; 5739.3432]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 0.7211 [0.3863; 2.2027]; tau = 0.8491 [0.6215; 1.4841]
## I^2 = 94.0% [91.8%; 95.6%]; H = 4.10 [3.50; 4.79]
##
## Test of heterogeneity:
## Q d.f. p-value
## 268.58 16 < 0.0001
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Log transformation
## - Normal approximation confidence interval for individual studies
## - Events per 100000 person-years##pooled rate = 3744 per 100,000 person-yearsPlots for reoffending
forest(forensic_patient_meta_result_reoffend, backtransf = TRUE,
xlab = "Incidence rate per 100,000 person-years",
refline = exp(forensic_patient_meta_result_violent$b),
addfit = TRUE,
xlim = c(-200, 35000))
funnel(forensic_patient_meta_result_reoffend)
funnel(forensic_patient_meta_result_reoffend,
xlim = c(-8, 8),
xlab = "Log Incidence Rate",
contour = c(0.90, 0.95, 0.99),
col.contour = c("gray80", "gray60", "gray40"),
backtransf = FALSE)
legend("topright", # position
legend = c("p < 0.01", "p < 0.05", "p < 0.10"),
fill = c("gray40", "gray60", "gray80"),
border = TRUE,
bty = "o",
title = "Significance regions",
cex = 0.9)
transformation for meta regression - reoffending
##<age>
#print(sort(forensic_patient_reoffend$m_age))
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(age_moderator = case_when(
!is.na(m_age) & m_age < 37.5 ~ 0,
!is.na(m_age) & m_age >= 37.5 ~ 1,
TRUE ~ NA_real_
))
#table(forensic_patient_reoffend$age_moderator)
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(ssd_percent = if_else(!is.na(ssd) & !is.na(sample), #checks for fields with no numbers
(ssd / sample) * 100,
NA_real_))
##<ssd>
#print(sort(forensic_patient_reoffend$ssd_percent))
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 70, 1, 0)))
#table(forensic_patient_reoffend$ssd_high)
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(ssd_percent = if_else(!is.na(ssd) & !is.na(sample), #checks for fields with no numbers
(ssd / sample) * 100,
NA_real_))
print(sort(forensic_patient_reoffend$ssd_percent))## [1] 29.54545 32.60870 45.71429 47.25275 53.19149 70.62500 71.34831 79.41176
## [9] 82.90155forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 70, 1, 0)))
#table(forensic_patient_reoffend$ssd_high, useNA = "ifany")
##<gender>
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(men_percent = if_else(!is.na(men) & !is.na(sample),
(men / sample) * 100,
NA_real_))
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(men_moderator = if_else(is.na(men_percent), NA_real_,
if_else(men_percent >= 90, 1, 0)))
#print(sort(forensic_patient_reoffend$men_moderator)) Meta-regression - reoffending
meta_reoffend <- forensic_patient_reoffend %>%
filter(!is.na(reoffending), !is.na(total_year_complete),
reoffending > 0, total_year_complete > 0) %>%
mutate(
rate = reoffending / total_year_complete * 100000, # per 100,000 person-years
yi = log(rate),
vi = 1 / reoffending # Poisson approximation for log incidence rate variance
)
# Step 3: Run meta-regression
meta_reoffend <- rma(yi = yi, vi = vi,
mods = ~ secure_level + ssd_high + m_follow,
data = meta_reoffend,
method = "REML")## Warning: 9 studies with NAs omitted from model fitting.## note - removed m_age, publication year, & gender because there weren't enough studies
summary(meta_reoffend)##
## Mixed-Effects Model (k = 8; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -5.8199 11.6398 21.6398 18.5713 81.6398
##
## tau^2 (estimated amount of residual heterogeneity): 0.7437 (SE = 0.6218)
## tau (square root of estimated tau^2 value): 0.8624
## I^2 (residual heterogeneity / unaccounted variability): 87.93%
## H^2 (unaccounted variability / sampling variability): 8.28
## R^2 (amount of heterogeneity accounted for): 7.80%
##
## Test for Residual Heterogeneity:
## QE(df = 4) = 18.2275, p-val = 0.0011
##
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 4.3643, p-val = 0.2247
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 8.8210 0.6752 13.0641 <.0001 7.4976 10.1444 ***
## secure_level 0.0942 0.3198 0.2946 0.7683 -0.5326 0.7210
## ssd_high -0.8269 0.8653 -0.9555 0.3393 -2.5228 0.8691
## m_follow -0.0891 0.1047 -0.8515 0.3945 -0.2942 0.1160
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1reconviction
#load excel file
forensic_patient_reconviction <- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 4)
View(forensic_patient_reconviction)
# <name variables to r friendly language>
forensic_patient_reconviction <- forensic_patient_reconviction %>%
rename(secure_level = `security level(0=unspecified or mixed/1=low/2=medium/3=high)`)
forensic_patient_reconviction <- forensic_patient_reconviction %>%
rename(after_2000 = `after_2000 (0=no/1=yes)`)random effect model for reconviction
forensic_patient_meta_result_reconvict <- metarate(event = reconviction,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_reconviction,
sm = "IRLN", # Log incidence rate
method = "Inverse", # DerSimonian-Laird (default) for random effects
common = FALSE, # Use random-effects model
irscale = 100000, # per 100,000 person-years
random = TRUE)
summary(forensic_patient_meta_result_reconvict)## events 95%-CI %W(random)
## Akande_2007 0.0000 [ 26.3253; 6728.7490] 1.2
## Halstead_2001 840.3361 [ 118.3727; 5965.6062] 2.0
## Cope_1993 2156.3342 [ 809.3107; 5745.3547] 3.7
## Simpson_2018 8333.3333 [3468.5651; 20021.0872] 3.9
## Nilsson_2011 3019.3237 [1256.7265; 7254.0171] 3.9
## Nagata_219 579.1506 [ 260.1896; 1289.1192] 4.1
## Simpson_2006 6944.4444 [3119.8659; 15457.4939] 4.1
## Blattner_2009 4186.2899 [2093.5537; 8370.9451] 4.3
## Dolan_2004 2795.0311 [1454.2959; 5371.8080] 4.4
## Haroon_2023 1621.3295 [ 843.6017; 3116.0551] 4.4
## Seifert_2005 693.0007 [ 360.5785; 1331.8874] 4.4
## Falla_2000 4844.2907 [2869.0427; 8179.4363] 4.7
## Davison_1999 6562.5000 [4278.7993; 10065.0682] 4.8
## Hogan_2019 3420.5830 [2273.0675; 5147.4000] 4.9
## Tabita_2012 2901.3540 [1944.6880; 4328.6403] 4.9
## Nicholls_2004 9784.0756 [6799.1659; 14079.3938] 4.9
## Thomson_2023 1172.7373 [ 837.9554; 1641.2722] 5.0
## Bailey_1992 6101.1905 [4492.4085; 8286.0954] 5.0
## Yoshikawa_2007 984.6247 [ 750.2920; 1292.1447] 5.0
## Friendship_1999 5336.3827 [4106.7700; 6934.1552] 5.1
## Ojansuu_2023 2015.0522 [1625.0053; 2498.7212] 5.1
## Salem_2015 6428.8559 [5346.3746; 7730.5073] 5.1
## Coid _2007 6255.4937 [5689.4340; 6877.8725] 5.2
##
## Number of studies: k = 23
## Number of events: e = 979
##
## events 95%-CI
## Random effects model 3024.8697 [2138.3665; 4278.8907]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 0.6009 [0.3219; 1.4022]; tau = 0.7751 [0.5673; 1.1842]
## I^2 = 94.6% [93.0%; 95.8%]; H = 4.30 [3.78; 4.90]
##
## Test of heterogeneity:
## Q d.f. p-value
## 407.60 22 < 0.0001
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Log transformation
## - Normal approximation confidence interval for individual studies
## - Continuity correction of 0.5 in studies with zero cell frequencies
## - Events per 100000 person-years##pooled rate = 3024 per 100,000 person-yearsPlots for reconviction
forest(forensic_patient_meta_result_reconvict, backtransf = TRUE,
xlab = "Incidence rate per 100,000 person-years",
refline = exp(forensic_patient_meta_result_violent$b),
addfit = TRUE,
xlim = c(-200, 35000))
funnel(forensic_patient_meta_result_reconvict,
xlim = c(-8, 8),
xlab = "Log Incidence Rate",
contour = c(0.90, 0.95, 0.99),
col.contour = c("gray80", "gray60", "gray40"),
backtransf = FALSE)
legend("topright", # position
legend = c("p < 0.01", "p < 0.05", "p < 0.10"),
fill = c("gray40", "gray60", "gray80"),
border = TRUE,
bty = "o",
title = "Significance regions",
cex = 0.9)
## transformation for meta regression - reconviction
## <Age>
#print(sort(forensic_patient_reconviction$m_age))
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(age_moderator = case_when(
!is.na(m_age) & m_age < 36.9 ~ 0,
!is.na(m_age) & m_age >= 37~ 1,
TRUE ~ NA_real_
))
#table(forensic_patient_reconviction$age_moderator)
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(ssd_percent = if_else(!is.na(ssd) & !is.na(sample), #checks for fields with no numbers
(ssd / sample) * 100,
NA_real_))
## <ssd>
#print(sort(forensic_patient_reconviction$ssd_percent))
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 70, 1, 0)))
#table(forensic_patient_reconviction$ssd_high)
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(ssd_percent = if_else(!is.na(ssd) & !is.na(sample), #checks for fields with no numbers
(ssd / sample) * 100,
NA_real_))
#print(sort(forensic_patient_reconviction$ssd_percent))
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 70, 1, 0)))
table(forensic_patient_reconviction$ssd_high, useNA = "ifany")##
## 0 1 <NA>
## 6 7 10forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(men_percent = if_else(!is.na(men) & !is.na(sample),
(men / sample) * 100,
NA_real_))
## <gender>
#print(sort(forensic_patient_reconviction$men_percent))
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(men_moderator = if_else(is.na(men_percent), NA_real_,
if_else(men_percent >= 84.9, 1, 0)))
#table(forensic_patient_reconviction$men_moderator, useNA = "ifany")Meta-regression - reconviction
meta_reconvict <- forensic_patient_reconviction %>%
filter(!is.na(reconviction), !is.na(total_year_complete),
reconviction > 0, total_year_complete > 0) %>%
mutate(
rate = reconviction / total_year_complete * 100000, # per 100,000 person-years
yi = log(rate),
vi = 1 / reconviction # Poisson approximation for log incidence rate variance
)
# Step 3: Run meta-regression
meta_reconvict <- rma(yi = yi, vi = vi,
mods = ~ secure_level + ssd_high + m_follow + after_2000 + m_age + men_moderator,
data = meta_reconvict,
method = "REML")## Warning: 14 studies with NAs omitted from model fitting.## Warning: Redundant predictors dropped from the model.summary(meta_reconvict)##
## Mixed-Effects Model (k = 8; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -2.4529 4.9058 18.9058 9.7578 130.9058
##
## tau^2 (estimated amount of residual heterogeneity): 0.2170 (SE = 0.3613)
## tau (square root of estimated tau^2 value): 0.4658
## I^2 (residual heterogeneity / unaccounted variability): 66.13%
## H^2 (unaccounted variability / sampling variability): 2.95
## R^2 (amount of heterogeneity accounted for): 52.71%
##
## Test for Residual Heterogeneity:
## QE(df = 2) = 5.0303, p-val = 0.0809
##
## Test of Moderators (coefficients 2:6):
## QM(df = 5) = 9.9219, p-val = 0.0775
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 14.2983 5.6328 2.5384 0.0111 3.2582 25.3384 *
## secure_level -0.4371 0.4563 -0.9578 0.3382 -1.3314 0.4573
## ssd_high 0.9212 1.6581 0.5556 0.5785 -2.3286 4.1710
## m_follow -0.1006 0.0499 -2.0168 0.0437 -0.1983 -0.0028 *
## m_age -0.1607 0.1600 -1.0044 0.3152 -0.4742 0.1529
## men_moderator 1.0174 0.7781 1.3075 0.1911 -0.5077 2.5425
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##colnames(model.matrix(meta_reg_reoffend))readmissions
#load excel file
forensic_patient_rehosp<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 5)
#View(forensic_patient_rehosp)
#glimpse(forensic_patient_rehosp)random effect model for redamission
forensic_patient_meta_result_rehosp <- metarate(event = readmission,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_rehosp,
sm = "IRLN", # Log incidence rate
method = "Inverse", # DerSimonian-Laird (default) for random effects
common = FALSE, # Use random-effects model
irscale = 100000, # per 100,000 person-years
random = TRUE)
summary(forensic_patient_meta_result_rehosp)## events 95%-CI %W(random)
## Blattner_2009 1046.5725 [ 261.7453; 4184.6553] 2.7
## Akande_2007 4208.7542 [ 1751.8006; 10111.6602] 3.7
## Halstead_2001 5882.3529 [ 2804.3172; 12338.8596] 3.9
## Simpson_2006 8101.8519 [ 3862.4276; 16994.4941] 3.9
## deVogel_2019 977.8357 [ 508.7824; 1879.3157] 4.1
## Falla_2000 3114.1869 [ 1620.3574; 5985.1979] 4.1
## Simpson_2018 25000.0000 [15071.6471; 41468.5932] 4.3
## Haroon_2023 3062.5113 [ 1903.8424; 4926.3401] 4.3
## Cope_1993 10781.6712 [ 6955.8728; 16711.6962] 4.4
## Penney_2018 27586.2069 [18490.1827; 41156.9115] 4.4
## Dolan_2004 10248.4472 [ 7285.8948; 14415.6171] 4.5
## Rossetto_2021 7070.7071 [ 5225.4007; 9567.6680] 4.6
## Tellefsen_1992 14000.0000 [10346.2933; 18943.9826] 4.6
## Jewell_2018 19371.5024 [14463.5090; 25944.9560] 4.6
## Hill_2024 4186.4139 [ 3198.3106; 5479.7871] 4.6
## Baxter_1999 22792.0228 [17540.2702; 29616.2087] 4.6
## Argo_2023 10819.3277 [ 8919.2615; 13124.1642] 4.7
## Lyons_2022* 7284.1444 [ 6072.2056; 8737.9716] 4.7
## Marshall_2014 18258.4270 [15867.4650; 21009.6669] 4.7
## Nagata_219 19787.6448 [17256.1382; 22690.5279] 4.7
## Salem_2015 16612.6188 [14812.3954; 18631.6322] 4.7
## Fovet_2022 12326.6857 [11413.8915; 13312.4780] 4.7
## Fazel_2016 4406.8665 [ 4279.5812; 4537.9376] 4.7
##
## Number of studies: k = 23
## Number of events: e = 6418
##
## events 95%-CI
## Random effects model 8567.5748 [6022.6965; 12187.7863]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 0.6853 [0.4011; 1.6003]; tau = 0.8279 [0.6334; 1.2650]
## I^2 = 98.9% [98.7%; 99.1%]; H = 9.50 [8.79; 10.28]
##
## Test of heterogeneity:
## Q d.f. p-value
## 1987.08 22 0
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Log transformation
## - Normal approximation confidence interval for individual studies
## - Events per 100000 person-years##pooled rate = 8567 per 100,000 person-yearsPlots
forest(forensic_patient_meta_result_rehosp, backtransf = TRUE,
addfit = TRUE)
funnel(forensic_patient_meta_result_rehosp,
xlim = c(-7, 7),
xlab = "Log Incidence Rate",
contour = c(0.90, 0.95, 0.99),
col.contour = c("gray80", "gray60", "gray40"),
backtransf = FALSE)
legend("topright", # position
legend = c("p < 0.01", "p < 0.05", "p < 0.10"),
fill = c("gray40", "gray60", "gray80"),
border = TRUE,
bty = "o",
title = "Significance regions",
cex = 0.9)
All-cause mortality
#load excel file
forensic_patient_mortality<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 6)
#View(forensic_patient_mortality)
#glimpse(forensic_patient_mortality)random effect model for all-cause mortality
forensic_patient_meta_result_mortality <- metarate(event = all_cause_mortality,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_mortality,
sm = "IRLN", # Log incidence rate
method = "Inverse", # DerSimonian-Laird (default) for random effects
common = FALSE, # Use random-effects model
irscale = 100000, # per 100,000 person-years
random = TRUE)
summary(forensic_patient_meta_result_mortality)## events 95%-CI %W(random)
## Nicholls_2004 1686.9096 [ 702.1387; 4052.8517] 7.2
## Argo_2023 735.2941 [ 350.5397; 1542.3574] 8.4
## Russo_1994 1538.4615 [ 733.4368; 3227.0863] 8.4
## deVogel_2019 1521.0778 [ 900.8620; 2568.2932] 10.7
## Takeda_2019 812.9304 [ 505.3667; 1307.6757] 11.2
## Tabita_2012 2417.7950 [1559.8578; 3747.6060] 11.6
## Thomson_2023 3173.2892 [2586.8173; 3892.7234] 13.7
## Ojansuu_2019 2780.8327 [2505.7285; 3086.1406] 14.3
## Fazel_2016 1920.6133 [1837.2112; 2007.8016] 14.4
##
## Number of studies: k = 9
## Number of events: e = 2465
##
## events 95%-CI
## Random effects model 1769.9378 [1268.4254; 2469.7392]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 0.1995 [0.0610; 0.8675]; tau = 0.4467 [0.2470; 0.9314]
## I^2 = 90.4% [84.0%; 94.2%]; H = 3.23 [2.50; 4.16]
##
## Test of heterogeneity:
## Q d.f. p-value
## 83.21 8 < 0.0001
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Log transformation
## - Normal approximation confidence interval for individual studies
## - Events per 100000 person-years##pooled rate = 1770 per 100,000 person-years##Plot for all-cause mortality
forest(forensic_patient_meta_result_mortality, backtransf = TRUE,
xlab = "Incidence rate per person-year")
funnel(forensic_patient_meta_result_mortality,
xlim = c(-7, 7),
xlab = "Log Incidence Rate",
contour = c(0.90, 0.95, 0.99),
col.contour = c("gray80", "gray60", "gray40"),
backtransf = FALSE)
legend("topright", # position
legend = c("p < 0.01", "p < 0.05", "p < 0.10"),
fill = c("gray40", "gray60", "gray80"),
border = TRUE,
bty = "o",
title = "Significance regions",
cex = 0.9)
suicide
#load excel file
forensic_patient_suicide<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 7)
#View(forensic_patient_suicide)
#glimpse(forensic_patient_suicide)random effect model for suicide
forensic_patient_meta_result_suicide <- metarate(event = suicide,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_suicide,
sm = "IRLN", # Log incidence rate
method = "Inverse", # DerSimonian-Laird (default) for random effects
common = FALSE, # Use random-effects model
irscale = 100000, # per 100,000 person-years
random = TRUE)
summary(forensic_patient_meta_result_suicide)## events 95%-CI %W(random)
## Argo_2023 210.0840 [ 52.5415; 840.0080] 4.1
## Russo_1994 439.5604 [109.9330; 1757.5552] 4.1
## deVogel_2019 325.9452 [105.1243; 1010.6158] 5.9
## Tabita_2012 725.3385 [325.8661; 1614.5158] 10.2
## Takeda_2019 478.1943 [257.2947; 888.7467] 14.3
## Ojansuu_2019 259.2302 [184.2936; 364.6370] 25.2
## Fazel_2016 436.5478 [397.7315; 479.1524] 36.2
##
## Number of studies: k = 7
## Number of events: e = 499
##
## events 95%-CI
## Random effects model 389.6390 [288.9856; 525.3498]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 0.0620 [0.0000; 0.6035]; tau = 0.2491 [0.0000; 0.7768]
## I^2 = 47.5% [0.0%; 77.8%]; H = 1.38 [1.00; 2.12]
##
## Test of heterogeneity:
## Q d.f. p-value
## 11.43 6 0.0759
##
## Details of meta-analysis methods:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tau
## - Calculation of I^2 based on Q
## - Log transformation
## - Normal approximation confidence interval for individual studies
## - Events per 100000 person-years##pooled rate = 390 per 100,000 person-yearssuicide outcome Plots
forest(forensic_patient_meta_result_suicide, backtransf = TRUE,
xlab = "Incidence rate per person-year",
xlim = c(-200, 2000))
funnel(forensic_patient_meta_result_suicide,
xlim = c(-8, 8),
xlab = "Log Incidence Rate",
contour = c(0.90, 0.95, 0.99),
col.contour = c("gray80", "gray60", "gray40"),
backtransf = FALSE)
legend("topright", # position
legend = c("p < 0.01", "p < 0.05", "p < 0.10"),
fill = c("gray40", "gray60", "gray80"),
border = TRUE,
bty = "o",
title = "Significance regions",
cex = 0.9)