fp_outcome_meta_analysis
James McLauchlan
2025-07-04
Demographic variables
#load excel file
forensic_patient_all<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 1)
# total number of studies
num_studies <- nrow(forensic_patient_all)
cat("Number of studies:", num_studies, "\n")## Number of studies: 48# total n
total_sample <- sum(forensic_patient_all$sample[forensic_patient_all$study_name != "Lyons_2022"], na.rm = TRUE)
cat("Total n(excluding studies using the same dataset):", total_sample, "\n")## Total n(excluding studies using the same dataset): 18318# security level
forensic_patient_all$secure_level_label <- factor(
forensic_patient_all$secure_level,
levels = 0:3,
labels = c("unspecified", "low", "medium", "high")
)
# n of studies reported total time at-risk
total_reported <- sum(forensic_patient_all$reported_at_risk == 1, na.rm = TRUE)
total_reported## [1] 6total_kaplan<- sum(forensic_patient_all$Reported_KM == 1, na.rm = TRUE)
total_kaplan## [1] 15table(forensic_patient_all$secure_level_label)##
## unspecified low medium high
## 29 2 11 6Violent reoffending
#load excel file
forensic_patient_violent<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 2)
forensic_patient_rob<- read_excel("fp_outcome_RoB.xlsx", sheet = 1)## New names:
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`forensic_patient_rob <- forensic_patient_rob %>%
mutate(study_name = stri_replace_all_regex(study_name, "\\p{Zs}+", "_"))
forensic_patient_violent <- forensic_patient_violent %>%
left_join(forensic_patient_rob %>%
select(study_name, total_score) %>%
rename(rob_score = total_score),
by = "study_name")
#view(forensic_patient_violent)
#View(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## 19 STUDIES WITH A COMBINED TOTAL OF 3352 reoffending incidents
## pooled incidence rate = 2375 per 100,000 person-years [1399.5739; 4033.3148]
## I^2 = 95.7% - approx. 96% of variability in study result is not due to chance - likely true differences between studies
## H = 4.80 - Significant sample variability (>1 = hetereogenity)
## Q = 414. p < .0001 - Q-test assesses whether observed variability between studies is greater than expected by chance (significant = p < .0001) Violent reoffending plots
forensic_patient_meta_result_violent$studlab <- gsub("^(.*)_(\\d{4})$", "\\1 et al. (\\2)", forensic_patient_meta_result_violent$studlab)
# Option 1: original forest plot
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,
prediction = TRUE,
xlim = c(-200, 35000))
# Option 1.5: Original plot (+RoB scores)
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,
prediction = TRUE,
xlim = c(-200, 35000),
leftcols = c("studlab", "effect", "ci", "rob_score"),
leftlabs = c("Study", "Events per 100,000", "95% CI", "RoB score"),
rightcols = FALSE
)
# Option 2: Cochrane review manager 5
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,
layout = "RevMan5",
xlim = c(-200, 35000),
leftcols = c("studlab", "effect", "ci"), # only include desired columns
leftlabs = c("Study", "Events per 100,000"),
prediction = TRUE,
rightcols = FALSE # removes any right-hand columns
)
#3: Journal of the American Medical Association guidelines
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,
layout = "JAMA",
xlim = c(-200, 35000),
prediction = TRUE,
leftcols = c("studlab", "effect", "ci"), # only include desired columns
leftlabs = c("Study", "Events per 100,000", "95% CI"),
rightcols = FALSE # removes any right-hand columns
)
funnel(forensic_patient_meta_result_violent,
xlim = c(-10, 10))
title("Funnel Plot (Violent Reoffending Studies)")
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)
graphics::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)
title("Contour-enhanced Funnel Plot (Violent Reoffending Studies)")
## Plot asymmetry suggest large variation in the true incidence rates across studies, even after accounting for within-study error
## There are no natural zero to compare against in prevalence studies, so the plot isn't anchored to a meaningful test of effects.
## colour-contour method might also not be conceptually meaningful because they're base don p-values testing whether the incidenceEgger’s test
##
metabias(forensic_patient_meta_result_violent, method.bias = "linreg")## Linear regression test of funnel plot asymmetry
##
## Test result: t = -0.02, df = 17, p-value = 0.9814
## Bias estimate: -0.0326 (SE = 1.3820)
##
## Details:
## - multiplicative residual heterogeneity variance (tau^2 = 24.4026)
## - predictor: standard error
## - weight: inverse variance
## - reference: Egger et al. (1997), BMJ## Eggers' test indicates no presence of funnel plot asymmetry, p = .98
## Funnel test is developed for effect sizes that compare two groups (e.g., odds ratio or the mean differences) where the expected true effect is usually assumed to be zero in the absence of treatment effect. Incidence rates are never expected to be zero, so there is no natural null effect to be tested against, and the asymmetry might be because the rates are low and skewed, not bias. Trim-and-fill method
# Run Trim and Fill (no backtransf argument here)
forensic_patient_trimfill <- trimfill(forensic_patient_meta_result_violent)
# View summary — make sure backtransf is a *logical* (TRUE or FALSE, not a string)
summary(forensic_patient_trimfill)## events 95%-CI %W(random)
## Coid et al. (2007) 3618.5174 [ 3194.2582; 4099.1266] 5.3
## Bengtson et al. (2019) 2380.9524 [ 2046.8173; 2769.6337] 5.3
## Baxter et al. (1999) 16687.0167 [12286.9291; 22662.8251] 5.3
## Dean et al. (2021) 144.6655 [ 54.2955; 385.4478] 4.6
## Blattner et al. (2009) 3139.7174 [ 1410.5516; 6988.6314] 4.9
## Fazel et al. (2016) 2574.9424 [ 2478.0819; 2675.5888] 5.3
## Thomson et al. (2023) 413.9073 [ 235.0620; 728.8259] 5.1
## Simpson et al. (2006) 6944.4444 [ 3119.8659; 15457.4939] 4.9
## Dowset et al. (2005) 6808.5106 [ 3404.9201; 13614.3626] 5.0
## Hogan et al. (2019) 2082.0940 [ 1233.1251; 3515.5519] 5.1
## deVogel et al. (2019) 1521.0778 [ 900.8620; 2568.2932] 5.1
## Miraglia et al. (2011) 740.1925 [ 542.9475; 1009.0936] 5.3
## Sivak et al. (2023) 3993.0556 [ 3153.7897; 5055.6614] 5.3
## Dolan et al. (2004) 1863.3540 [ 837.1317; 4147.6008] 4.9
## Simpson et al. (2018) 16666.6667 [ 8967.5785; 30975.7843] 5.0
## Ojansuu et al. (2023) 1165.3314 [ 878.1912; 1546.3571] 5.3
## Edwards et al. (2002) 303.0303 [ 42.6859; 2151.2338] 3.3
## Davison et al. (1999) 5000.0000 [ 3063.1596; 8161.5076] 5.2
## Friendship et al. (1999) 2763.4839 [ 1920.4048; 3976.6841] 5.2
## Filled: Dean et al. (2021) 48646.3376 [18257.8392; 129613.7038] 4.6
##
## Number of studies: k = 20 (with 1 added studies)
##
## events 95%-CI
## Random effects model 2711.3952 [1522.7036; 4828.0334]
##
## Quantifying heterogeneity (with 95%-CIs):
## tau^2 = 1.6198 [0.9112; 3.9202]; tau = 1.2727 [0.9546; 1.9799]
## I^2 = 95.8% [94.5%; 96.7%]; H = 4.86 [4.27; 5.53]
##
## Test of heterogeneity:
## Q d.f. p-value
## 448.75 19 < 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
## - Trim-and-fill method to adjust for funnel plot asymmetry (L-estimator)
## - Log transformation
## - Events per 100000 person-years# Funnel plot
funnel(forensic_patient_trimfill,
xlab = "Log incidence rate",
col = "blue",
backtransf = FALSE)
title("Trim-and-fill method (Violent reoffending studies)")
## Added only one study representing a hypothetical missing finding was imputed by the algorithm to balance the funnel plot - minimal funnel plot asymmetry
## I^2 = 95.8% - Most variation in effect sizes is due to between-study heterogeneity, not chance
## Adjusted pooled incidence rate = 2711 per 100,000 person-years [1522.7036; 4828.0334]Rucker’s limit meta-analysis method
res_limit <- limitmeta(forensic_patient_meta_result_violent)
summary(res_limit)## Results for individual studies
## (left: original data; right: shrunken estimates)
##
## IRLN 95%-CI IRLN 95%-CI
## Coid et al. (2007) 0.0362 [0.0319; 0.0410] 0.0366 [0.0323; 0.0415]
## Bengtson et al. (2019) 0.0238 [0.0205; 0.0277] 0.0242 [0.0208; 0.0282]
## Baxter et al. (1999) 0.1669 [0.1229; 0.2266] 0.1764 [0.1299; 0.2396]
## Dean et al. (2021) 0.0014 [0.0005; 0.0039] 0.0036 [0.0013; 0.0096]
## Blattner et al. (2009) 0.0314 [0.0141; 0.0699] 0.0493 [0.0221; 0.1097]
## Fazel et al. (2016) 0.0257 [0.0248; 0.0268] 0.0258 [0.0248; 0.0268]
## Thomson et al. (2023) 0.0041 [0.0024; 0.0073] 0.0056 [0.0032; 0.0098]
## Simpson et al. (2006) 0.0694 [0.0312; 0.1546] 0.1040 [0.0467; 0.2314]
## Dowset et al. (2005) 0.0681 [0.0340; 0.1361] 0.0929 [0.0464; 0.1857]
## Hogan et al. (2019) 0.0208 [0.0123; 0.0352] 0.0258 [0.0153; 0.0436]
## deVogel et al. (2019) 0.0152 [0.0090; 0.0257] 0.0190 [0.0113; 0.0321]
## Miraglia et al. (2011) 0.0074 [0.0054; 0.0101] 0.0081 [0.0059; 0.0110]
## Sivak et al. (2023) 0.0399 [0.0315; 0.0506] 0.0416 [0.0329; 0.0527]
## Dolan et al. (2004) 0.0186 [0.0084; 0.0415] 0.0302 [0.0136; 0.0671]
## Simpson et al. (2018) 0.1667 [0.0897; 0.3098] 0.2074 [0.1116; 0.3855]
## Ojansuu et al. (2023) 0.0117 [0.0088; 0.0155] 0.0125 [0.0094; 0.0166]
## Edwards et al. (2002) 0.0030 [0.0004; 0.0215] 0.0366 [0.0051; 0.2595]
## Davison et al. (1999) 0.0500 [0.0306; 0.0816] 0.0591 [0.0362; 0.0965]
## Friendship et al. (1999) 0.0276 [0.0192; 0.0398] 0.0306 [0.0213; 0.0440]
##
## Result of limit meta-analysis:
##
## Random effects model events 95%-CI z pval
## Adjusted estimate 3168.3916 [1691.4288; 5935.0444] -10.78 < 0.0001
## Unadjusted estimate 2375.9045 [1399.5739; 4033.3148] -- --
##
## Quantifying heterogeneity:
## tau^2 = 1.2817; I^2 = 95.7% [94.3%; 96.7%]; G^2 = 1.3%
##
## Test of heterogeneity:
## Q d.f. p-value
## 414.86 18 0
##
## Test of small-study effects:
## Q-Q' d.f. p-value
## 0.01 1 0.9071
##
## Test of residual heterogeneity beyond small-study effects:
## Q' d.f. p-value
## 414.84 17 0
##
## Adjustment method: expectation (beta0)# Q-Q = 0.01 p = .91 no evidence of small-study effects or publication bias
# the model is a better fit when accounting for heterogeneity
## I^2 = 95.7% - high between-study differences
## G^2 = 1.3% = only 1.3% of heterogeneity is explainde by small-study effects, likely not a major source of variability
## adjusted pooled incidence rate = 3168 per 100,000 person-years [[1691.4288; 5935.0444]
## The difference might reflect what the effect sizes would be if the studies had small standard error - i.e., larger studies (9 out of 18 studies included had n < 100)
# Funnel with curve
funnel.limitmeta(res_limit)
# Funnel with curve and shrunken study estimates
funnel.limitmeta(res_limit, shrunken = TRUE)
## funnel plot 1 grey curve indicates average effect size when standard error on y-axis is 0
## funnel plot 2 helps visualise how the studies would look if those biases were accounted fortransformation 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))
forensic_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")
#< publication < 15yrs>
forensic_patient_violent <- forensic_patient_violent %>%
mutate(published_year = as.numeric(str_extract(study_name, "\\d{4}")))
#str(forensic_patient_violent$published_year)
#print(sort(forensic_patient_violent$published_year))
forensic_patient_violent <- forensic_patient_violent %>%
mutate(published_after_2010 = ifelse(published_year >= 2010, 1, 0))
#print(sort(forensic_patient_violent$published_after_2010))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
)
res <- rma(yi = yi, vi = vi,
mods = ~ secure_level + published_after_2010 + age_moderator + ssd_percent + uk_studies + m_follow,
data = forensic_patient_violent,
method = "REML")## Warning: 9 studies with NAs omitted from model fitting.summary(res)##
## Mixed-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -3.0377 6.0754 22.0754 14.8643 166.0754
##
## tau^2 (estimated amount of residual heterogeneity): 0.4096 (SE = 0.3760)
## tau (square root of estimated tau^2 value): 0.6400
## I^2 (residual heterogeneity / unaccounted variability): 91.45%
## H^2 (unaccounted variability / sampling variability): 11.69
## R^2 (amount of heterogeneity accounted for): 59.49%
##
## Test for Residual Heterogeneity:
## QE(df = 3) = 40.2251, p-val < .0001
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 17.6213, p-val = 0.0073
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 11.6244 1.2251 9.4883 <.0001 9.2232 14.0256 ***
## secure_level 0.0148 0.2231 0.0664 0.9471 -0.4225 0.4521
## published_after_2010 0.5905 1.2616 0.4680 0.6398 -1.8823 3.0632
## age_moderator -0.7410 0.9183 -0.8070 0.4197 -2.5408 1.0588
## ssd_percent -0.0252 0.0223 -1.1274 0.2596 -0.0689 0.0186
## uk_studies -0.5033 0.9873 -0.5098 0.6102 -2.4384 1.4317
## m_follow -0.2030 0.0670 -3.0315 0.0024 -0.3342 -0.0718 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## A longer follow-up duration is significantly associated with a lower violent reoffending rates
# sensitivity analysis with moderators close to significant value: m follow-up & ssd percent
res_reduced <- rma(yi = yi, vi = vi,
mods = ~ ssd_percent + m_follow,
data = forensic_patient_violent,
method = "REML")## Warning: 9 studies with NAs omitted from model fitting.summary(res_reduced)##
## Mixed-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -6.7083 13.4167 21.4167 21.2003 41.4167
##
## tau^2 (estimated amount of residual heterogeneity): 0.3311 (SE = 0.2059)
## tau (square root of estimated tau^2 value): 0.5754
## I^2 (residual heterogeneity / unaccounted variability): 94.03%
## H^2 (unaccounted variability / sampling variability): 16.76
## R^2 (amount of heterogeneity accounted for): 67.25%
##
## Test for Residual Heterogeneity:
## QE(df = 7) = 96.8444, p-val < .0001
##
## Test of Moderators (coefficients 2:3):
## QM(df = 2) = 17.2816, p-val = 0.0002
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 10.9438 0.9752 11.2223 <.0001 9.0325 12.8551 ***
## ssd_percent -0.0282 0.0119 -2.3697 0.0178 -0.0515 -0.0049 *
## m_follow -0.1318 0.0335 -3.9341 <.0001 -0.1975 -0.0662 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## The model accounted for 67.25% of the between-study heterogeneity (R² = 67.25%) and explained significantly more variance than expected by chance (QM(2) = 17.28, p = .0002).
## ssd_percent, although non-significant in the full model, may represent a meaningful moderator when unrelated variables are excluded
## higher proportion of participants with schizophrenia-spectrum disorders was associated with significantly lower rates of violent reoffending (β = −0.0282, p = .0178), and longer follow-up periods were similarly associated with reduced rates (β = −0.1318, p < .0001).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,
prediction = TRUE,
xlim = c(-200, 25000))
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)
graphics::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))
#< publication < 15yrs>
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(published_year = as.numeric(str_extract(study_name, "\\d{4}")))
#print(sort(forensic_patient_reoffend$published_year))
forensic_patient_reoffend <- forensic_patient_reoffend %>%
mutate(published_after_2010 = ifelse(published_year >= 2010, 1, 0))
#print(sort(forensic_patient_reoffend$published_after_2010))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
)
meta_reoffend <- rma(yi = yi, vi = vi,
mods = ~ secure_level + ssd_percent + m_follow + published_after_2010,
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
## -4.9027 9.8053 21.8053 16.3970 105.8053
##
## tau^2 (estimated amount of residual heterogeneity): 1.2053 (SE = 1.1033)
## tau (square root of estimated tau^2 value): 1.0978
## I^2 (residual heterogeneity / unaccounted variability): 92.39%
## H^2 (unaccounted variability / sampling variability): 13.14
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 3) = 17.9188, p-val = 0.0005
##
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 2.9970, p-val = 0.5583
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 9.4181 1.3489 6.9822 <.0001 6.7744 12.0619 ***
## secure_level 0.0951 0.4385 0.2169 0.8283 -0.7643 0.9544
## ssd_percent -0.0125 0.0293 -0.4268 0.6695 -0.0700 0.0450
## m_follow -0.0919 0.1333 -0.6893 0.4906 -0.3531 0.1694
## published_after_2010 -0.4618 1.1908 -0.3878 0.6981 -2.7958 1.8721
##
## ---
## 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-years
##continuity correction applied because the Akande et al. study reported 0 reconviction
res_reconviction1 <- metarate(event = reconviction,
time = total_year_complete,
studlab = study_name,
data = forensic_patient_reconviction,
sm = "IRLN", # Log incidence rate
method = "Inverse", # Use DerSimonian-Laird
random = TRUE,
irscale = 100000) # per 100,000 person-years
summary(res_reconviction1) ## outcome is the same## events 95%-CI %W(common) %W(random)
## Akande_2007 0.0000 [ 26.3253; 6728.7490] 0.1 1.2
## Halstead_2001 840.3361 [ 118.3727; 5965.6062] 0.1 2.0
## Cope_1993 2156.3342 [ 809.3107; 5745.3547] 0.4 3.7
## Simpson_2018 8333.3333 [3468.5651; 20021.0872] 0.5 3.9
## Nilsson_2011 3019.3237 [1256.7265; 7254.0171] 0.5 3.9
## Nagata_219 579.1506 [ 260.1896; 1289.1192] 0.6 4.1
## Simpson_2006 6944.4444 [3119.8659; 15457.4939] 0.6 4.1
## Blattner_2009 4186.2899 [2093.5537; 8370.9451] 0.8 4.3
## Dolan_2004 2795.0311 [1454.2959; 5371.8080] 0.9 4.4
## Haroon_2023 1621.3295 [ 843.6017; 3116.0551] 0.9 4.4
## Seifert_2005 693.0007 [ 360.5785; 1331.8874] 0.9 4.4
## Falla_2000 4844.2907 [2869.0427; 8179.4363] 1.4 4.7
## Davison_1999 6562.5000 [4278.7993; 10065.0682] 2.1 4.8
## Hogan_2019 3420.5830 [2273.0675; 5147.4000] 2.3 4.9
## Tabita_2012 2901.3540 [1944.6880; 4328.6403] 2.5 4.9
## Nicholls_2004 9784.0756 [6799.1659; 14079.3938] 3.0 4.9
## Thomson_2023 1172.7373 [ 837.9554; 1641.2722] 3.5 5.0
## Bailey_1992 6101.1905 [4492.4085; 8286.0954] 4.2 5.0
## Yoshikawa_2007 984.6247 [ 750.2920; 1292.1447] 5.3 5.0
## Friendship_1999 5336.3827 [4106.7700; 6934.1552] 5.7 5.1
## Ojansuu_2023 2015.0522 [1625.0053; 2498.7212] 8.5 5.1
## Salem_2015 6428.8559 [5346.3746; 7730.5073] 11.5 5.1
## Coid _2007 6255.4937 [5689.4340; 6877.8725] 43.6 5.2
##
## Number of studies: k = 23
## Number of events: e = 979
##
## events 95%-CI
## Common effect model 4413.0525 [4145.1621; 4698.2559]
## 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-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,
prediction = TRUE,
xlim = c(-200, 25000))
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)
graphics::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")
forensic_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")
#< publication < 15yrs>
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(published_year = as.numeric(str_extract(study_name, "\\d{4}")))
#print(sort(forensic_patient_reconviction$published_year))
forensic_patient_reconviction <- forensic_patient_reconviction %>%
mutate(published_after_2010 = ifelse(published_year >= 2010, 1, 0))
#print(sort(forensic_patient_reconviction$published_after_2010))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 + published_after_2010 + m_age + men_moderator,
data = meta_reconvict,
method = "REML")## Warning: 14 studies with NAs omitted from model fitting.summary(meta_reconvict)##
## Mixed-Effects Model (k = 8; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -1.2114 2.4228 18.4228 2.4228 162.4228
##
## tau^2 (estimated amount of residual heterogeneity): 0.5172 (SE = 0.9338)
## tau (square root of estimated tau^2 value): 0.7192
## I^2 (residual heterogeneity / unaccounted variability): 78.34%
## H^2 (unaccounted variability / sampling variability): 4.62
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 1) = 4.6158, p-val = 0.0317
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 5.5248, p-val = 0.4785
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 4.2547 12.2568 0.3471 0.7285 -19.7682 28.2776
## secure_level 0.4077 1.0734 0.3798 0.7041 -1.6962 2.5116
## ssd_high -2.4194 3.9610 -0.6108 0.5413 -10.1827 5.3439
## m_follow -0.1366 0.0835 -1.6371 0.1016 -0.3002 0.0269
## published_after_2010 1.3484 1.3613 0.9906 0.3219 -1.3197 4.0165
## m_age 0.0928 0.3279 0.2830 0.7772 -0.5499 0.7355
## men_moderator 1.5861 1.1653 1.3611 0.1735 -0.6979 3.8701
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# sensitivity analysis - removed mean age and gender as moderators
##colnames(model.matrix(meta_reg_reoffend))Meta-regression - reconviction - sensitivity analysis
meta_reconvict2 <- 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_reconvict3 <- rma(yi = yi, vi = vi,
mods = ~ secure_level + ssd_high + m_follow + published_after_2010,
data = meta_reconvict2,
method = "REML")## Warning: 12 studies with NAs omitted from model fitting.summary(meta_reconvict3)##
## Mixed-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## -4.6750 9.3500 21.3500 19.0066 105.3500
##
## tau^2 (estimated amount of residual heterogeneity): 0.1262 (SE = 0.1332)
## tau (square root of estimated tau^2 value): 0.3552
## I^2 (residual heterogeneity / unaccounted variability): 68.11%
## H^2 (unaccounted variability / sampling variability): 3.14
## R^2 (amount of heterogeneity accounted for): 63.90%
##
## Test for Residual Heterogeneity:
## QE(df = 5) = 14.0905, p-val = 0.0150
##
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 12.7848, p-val = 0.0124
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 8.7370 0.3667 23.8277 <.0001 8.0184 9.4557 ***
## secure_level 0.1275 0.1755 0.7264 0.4676 -0.2165 0.4716
## ssd_high -0.2608 0.3064 -0.8512 0.3947 -0.8614 0.3397
## m_follow -0.1035 0.0431 -2.3994 0.0164 -0.1880 -0.0190 *
## published_after_2010 0.1117 0.3959 0.2821 0.7779 -0.6643 0.8877
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# sensitivity analysis - removed mean age and gender as moderators
##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)
forensic_patient_rehosp <- forensic_patient_rehosp %>%
rename(secure_level = `security level(0=unspecified or mixed/1=low/2=medium/3=high)`)
forensic_patient_rehosp <- forensic_patient_rehosp %>%
rename(after_2000 = `after_2000 (0=no/1=yes)`)
str(forensic_patient_rehosp)## tibble [23 × 33] (S3: tbl_df/tbl/data.frame)
## $ study_name : chr [1:23] "Blattner_2009" "Akande_2007" "Halstead_2001" "Simpson_2006" ...
## $ country : chr [1:23] "uk" "uk" "uk" "nz" ...
## $ year : num [1:23] 2009 2007 2001 2006 2019 ...
## $ secure_level : num [1:23] 3 1 2 0 0 2 0 0 2 0 ...
## $ after_2000 : num [1:23] 1 1 1 1 1 1 1 1 0 1 ...
## $ Source (0 =systematic review/1 = received from authors and citation lists: num [1:23] 0 0 0 0 0 1 0 0 1 0 ...
## $ Reported KM curve (0 = no, 1 = yes) : num [1:23] 0 0 0 0 0 0 0 0 0 1 ...
## $ reported time at risk : num [1:23] 0 0 0 0 0 0 0 0 0 0 ...
## $ Reported PY : num [1:23] 0 0 0 0 0 0 0 0 0 0 ...
## $ sample : num [1:23] 39 33 35 48 78 85 60 61 35 87 ...
## $ m_age : num [1:23] 36.4 35.8 29.3 37.5 38.6 ...
## $ men : num [1:23] NA 25 29 NA 0 NA 51 50 NA 73 ...
## $ women : num [1:23] NA 8 6 48 71 NA 9 11 NA 14 ...
## $ trans : num [1:23] 0 0 0 0 0 NA 0 0 NA 0 ...
## $ ssd : num [1:23] NA 29 16 NA NA NA 49 NA NA 72 ...
## $ m_follow : num [1:23] 4.9 3.6 3.4 1.8 11.8 3.4 1 9.1 5.3 1 ...
## $ median_follow : num [1:23] NA NA NA NA NA NA NA NA NA NA ...
## $ total_yr_at_risk : num [1:23] NA NA NA NA NA NA NA NA NA NA ...
## $ violent_reoffend : num [1:23] 6 NA NA 6 14 NA 10 NA NA NA ...
## $ reoffending : num [1:23] NA NA 11 9 NA NA NA NA NA NA ...
## $ reconviction : num [1:23] 8 0 1 6 NA 14 5 9 4 NA ...
## $ readmission : num [1:23] 2 5 7 7 9 9 15 17 20 24 ...
## $ all_cause_mortality : num [1:23] NA NA NA NA 14 NA NA NA NA NA ...
## $ suicide : num [1:23] NA NA NA NA 3 NA NA NA NA NA ...
## $ violent_outcome : num [1:23] NA NA NA NA NA NA NA NA NA NA ...
## $ reoffend_outcome : logi [1:23] NA NA NA NA NA NA ...
## $ reconvict_outcome : logi [1:23] NA NA NA NA NA NA ...
## $ readmit_outcome : num [1:23] NA NA NA NA NA NA NA NA NA NA ...
## $ allcause_outcome : num [1:23] NA NA NA NA NA NA NA NA NA NA ...
## $ suicide_outcome : logi [1:23] NA NA NA NA NA NA ...
## $ mean_estimate : num [1:23] 191.1 118.8 119 86.4 920.4 ...
## $ median_estimate : num [1:23] NA NA NA NA NA NA NA NA NA NA ...
## $ total_year_complete : num [1:23] 191.1 118.8 119 86.4 920.4 ...Random effect model for redamissions
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_2019 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 for redamission
forest(forensic_patient_meta_result_rehosp, backtransf = TRUE,
xlab = "Incidence rate per 100,000 person-years",
refline = exp(forensic_patient_meta_result_rehosp$b),
addfit = TRUE,
prediction = TRUE,
xlim = c(-200, 55000))
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)
# Now add the legend
graphics::legend(
"topright",
legend = c("p < 0.10", "p < 0.05", "p < 0.01"),
col = c("gray80", "gray60", "gray40"),
bty = "o",
title = "Significance regions",
cex = 0.9
)
Transformation for meta regression - radmission
## <Age>
print(sort(forensic_patient_rehosp$m_age))## [1] 29.30 31.30 33.60 35.00 35.80 36.30 36.40 36.44 36.60 37.50 38.20 38.60
## [13] 39.62 40.00 44.00 44.50forensic_patient_rehosp <- forensic_patient_rehosp %>%
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_rehosp$age_moderator)
## <ssd>
forensic_patient_rehosp <- forensic_patient_rehosp %>%
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_rehosp$ssd_percent))
forensic_patient_rehosp <- forensic_patient_rehosp %>%
mutate(ssd_high = if_else(is.na(ssd_percent), NA_real_,
if_else(ssd_percent >= 70, 1, 0)))
#table(forensic_patient_rehosp$ssd_high)
#table(forensic_patient_rehosp$ssd_high, useNA = "ifany")
## <gender>
forensic_patient_rehosp <- forensic_patient_rehosp %>%
mutate(men_percent = if_else(!is.na(men) & !is.na(sample),
(men / sample) * 100,
NA_real_))
#print(sort(forensic_patient_rehosp$men_percent))
#< publication < 15yrs>
forensic_patient_rehosp <- forensic_patient_rehosp %>%
mutate(published_year = as.numeric(str_extract(study_name, "\\d{4}")))
#print(sort(forensic_patient_rehosp$published_year))
forensic_patient_rehosp <- forensic_patient_rehosp %>%
mutate(published_after_2010 = ifelse(published_year >= 2010, 1, 0))
#print(sort(forensic_patient_rehosp$published_after_2010))Meta-regression - readmission
meta_rehosp <- forensic_patient_rehosp %>%
filter(!is.na(readmission), !is.na(total_year_complete),
readmission > 0, total_year_complete > 0) %>%
mutate(
rate = readmission / total_year_complete * 100000, # per 100,000 person-years
yi = log(rate),
vi = 1 / readmission # Poisson approximation for log incidence rate variance
)
# Step 3: Run meta-regression
meta_rehosp <- rma(yi = yi, vi = vi,
mods = ~ secure_level + ssd_percent + m_follow + published_after_2010 + m_age + men_percent,
data = meta_rehosp,
method = "REML")## Warning: 13 studies with NAs omitted from model fitting.summary(meta_rehosp)##
## Mixed-Effects Model (k = 10; tau^2 estimator: REML)
##
## logLik deviance AIC BIC AICc
## 1.8404 -3.6808 12.3192 5.1081 156.3192
##
## tau^2 (estimated amount of residual heterogeneity): 0 (SE = 0.0113)
## tau (square root of estimated tau^2 value): 0
## I^2 (residual heterogeneity / unaccounted variability): 0.00%
## H^2 (unaccounted variability / sampling variability): 1.00
## R^2 (amount of heterogeneity accounted for): 100.00%
##
## Test for Residual Heterogeneity:
## QE(df = 3) = 2.2037, p-val = 0.5312
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 1022.9444, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 6.1797 1.0443 5.9177 <.0001 4.1330 8.2265 ***
## secure_level 0.7023 0.5491 1.2790 0.2009 -0.3739 1.7785
## ssd_percent 0.0061 0.0090 0.6714 0.5020 -0.0116 0.0238
## m_follow -0.0913 0.0248 -3.6789 0.0002 -0.1400 -0.0427 ***
## published_after_2010 2.0835 1.0410 2.0015 0.0453 0.0432 4.1237 *
## m_age 0.0240 0.0260 0.9219 0.3566 -0.0270 0.0751
## men_percent 0.0053 0.0024 2.2076 0.0273 0.0006 0.0099 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##colnames(model.matrix(meta_reg_reoffend))
## 13 studies were removed because of missing value in one or more variable
# 11 studies missing ssd percent - remove
##meta-regression: readmission (version 2)
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-yearsPlot for all-cause mortality
forest(forensic_patient_meta_result_mortality, backtransf = TRUE,
xlab = "Incidence rate per 100,000 person-years",
refline = exp(forensic_patient_meta_result_mortality$b),
addfit = TRUE,
prediction = TRUE,
xlim = c(-200, 7000))
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)
graphics::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)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 100,000 person-years",
refline = exp(forensic_patient_meta_result_suicide$b),
addfit = TRUE,
prediction = TRUE,
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)
graphics::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)