fp_outcome_meta_analysis
James McLauchlan
2025-07-04
Violent reoffending
#load excel file
forensic_patient_violent<- read_excel("fp_data_estimated_py_all_outcomes.xlsx", sheet = 2)
#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## pooled rate = 2375 per 100,000 person-yearsPrediction interval for violent reoffending
forensic_patient_violent$yi <- log(forensic_patient_violent$violent_reoffend / forensic_patient_violent$total_year_complete)
forensic_patient_violent$vi <- 1 / forensic_patient_violent$violent_reoffend
res_violent <- rma(yi = yi, vi = vi, data = forensic_patient_violent, method = "REML")
predict(res_violent , transf = exp) # Back-transforms to incidence rate scale##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.0238 0.0140 0.0403 0.0024 0.2325# Extract parameters
mu <- res_violent$b[1] # pooled log incidence rate
se <- res_violent$se[1] # standard error of pooled effect
tau2 <- res_violent$tau2 # between-study variance (tau squared)
# 95% CI of pooled effect (log scale)
ci_lb_log <- mu - 1.96 * se
ci_ub_log <- mu + 1.96 * se
# 95% Prediction interval (log scale)
pi_lb_log <- mu - 1.96 * sqrt(se^2 + tau2)
pi_ub_log <- mu + 1.96 * sqrt(se^2 + tau2)
# Back-transform to incidence rate scale
ci_lb <- exp(ci_lb_log)
ci_ub <- exp(ci_ub_log)
pi_lb <- exp(pi_lb_log)
pi_ub <- exp(pi_ub_log)
# Scale to per 100,000 person-years
ci_lb_100k <- ci_lb * 100000
ci_ub_100k <- ci_ub * 100000
pi_lb_100k <- pi_lb * 100000
pi_ub_100k <- pi_ub * 100000
pooled_rate_100k <- exp(mu) * 100000
# Print results
cat("Pooled incidence rate per 100,000 person-years:\n")## Pooled incidence rate per 100,000 person-years:cat(sprintf("%.2f [%.2f, %.2f]\n", pooled_rate_100k, ci_lb_100k, ci_ub_100k))## 2375.90 [1399.56, 4033.35]cat("\nPrediction interval per 100,000 person-years:\n")##
## Prediction interval per 100,000 person-years:cat(sprintf("[%.2f, %.2f]\n", pi_lb_100k, pi_ub_100k))## [242.74, 23255.34]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,
xlim = c(-200, 35000))
# 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", "95% CI"),
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),
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)
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 errorEgger’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 epresence of funnel plot asymmetryTrim-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)")
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
# Funnel with curve
funnel.limitmeta(res_limit)
# Funnel with curve and shrunken study estimates
funnel.limitmeta(res_limit, shrunken = TRUE)
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))
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
)
# Step 3: Run meta-regression
res <- rma(yi = yi, vi = vi,
mods = ~ secure_level + published_after_2010 + age_moderator + ssd_high + 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
## -2.6638 5.3276 21.3276 14.1165 165.3276
##
## tau^2 (estimated amount of residual heterogeneity): 0.2948 (SE = 0.2702)
## tau (square root of estimated tau^2 value): 0.5430
## I^2 (residual heterogeneity / unaccounted variability): 91.28%
## H^2 (unaccounted variability / sampling variability): 11.46
## R^2 (amount of heterogeneity accounted for): 70.84%
##
## Test for Residual Heterogeneity:
## QE(df = 3) = 23.1950, p-val < .0001
##
## Test of Moderators (coefficients 2:7):
## QM(df = 6) = 23.6744, p-val = 0.0006
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 9.4969 1.3346 7.1158 <.0001 6.8811 12.1127 ***
## secure_level 0.0606 0.1991 0.3042 0.7610 -0.3297 0.4509
## published_after_2010 1.5973 1.4481 1.1031 0.2700 -1.2408 4.4354
## age_moderator -0.2910 0.9402 -0.3095 0.7570 -2.1337 1.5518
## ssd_high -1.4108 0.9035 -1.5615 0.1184 -3.1816 0.3600
## uk_studies 0.1379 1.0516 0.1311 0.8957 -1.9233 2.1990
## m_follow -0.2072 0.0586 -3.5357 0.0004 -0.3220 -0.0923 ***
##
## ---
## 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-yearsPrediction interval for reoffending
forensic_patient_reoffend$yi <- log(forensic_patient_reoffend$reoffending / forensic_patient_reoffend$total_year_complete)
forensic_patient_reoffend$vi <- 1 / forensic_patient_reoffend$reoffending
res_reoffend <- rma(yi = yi, vi = vi, data = forensic_patient_reoffend, method = "REML")
predict(res_reoffend , transf = exp) # Back-transforms to incidence rate scale##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.0374 0.0244 0.0574 0.0067 0.2087# Extract parameters
mu <- res_reoffend$b[1] # pooled log incidence rate
se <- res_reoffend$se[1] # standard error of pooled effect
tau2 <- res_reoffend$tau2 # between-study variance (tau squared)
# 95% CI of pooled effect (log scale)
ci_lb_log <- mu - 1.96 * se
ci_ub_log <- mu + 1.96 * se
# 95% Prediction interval (log scale)
pi_lb_log <- mu - 1.96 * sqrt(se^2 + tau2)
pi_ub_log <- mu + 1.96 * sqrt(se^2 + tau2)
# Back-transform to incidence rate scale
ci_lb <- exp(ci_lb_log)
ci_ub <- exp(ci_ub_log)
pi_lb <- exp(pi_lb_log)
pi_ub <- exp(pi_ub_log)
# Scale to per 100,000 person-years
ci_lb_100k <- ci_lb * 100000
ci_ub_100k <- ci_ub * 100000
pi_lb_100k <- pi_lb * 100000
pi_ub_100k <- pi_ub * 100000
pooled_rate_100k <- exp(mu) * 100000
# Print results
cat("Pooled incidence rate per 100,000 person-years:\n")## Pooled incidence rate per 100,000 person-years:cat(sprintf("%.2f [%.2f, %.2f]\n", pooled_rate_100k, ci_lb_100k, ci_ub_100k))## 3744.74 [2443.30, 5739.39]cat("\nPrediction interval per 100,000 person-years:\n")##
## Prediction interval per 100,000 person-years:cat(sprintf("[%.2f, %.2f]\n", pi_lb_100k, pi_ub_100k))## [671.74, 20875.65]## Borenstein's 95% prediction interval 671 - 20,870 per 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, 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)
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_high + 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.8602 9.7205 21.7205 16.3122 105.7205
##
## tau^2 (estimated amount of residual heterogeneity): 1.1670 (SE = 1.0727)
## tau (square root of estimated tau^2 value): 1.0803
## I^2 (residual heterogeneity / unaccounted variability): 92.01%
## H^2 (unaccounted variability / sampling variability): 12.52
## R^2 (amount of heterogeneity accounted for): 0.00%
##
## Test for Residual Heterogeneity:
## QE(df = 3) = 18.1414, p-val = 0.0004
##
## Test of Moderators (coefficients 2:5):
## QM(df = 4) = 3.1606, p-val = 0.5313
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 8.8866 0.9537 9.3183 <.0001 7.0174 10.7557 ***
## secure_level 0.0679 0.4277 0.1589 0.8738 -0.7702 0.9061
## ssd_high -0.7554 1.4474 -0.5219 0.6017 -3.5923 2.0815
## m_follow -0.0878 0.1306 -0.6726 0.5012 -0.3437 0.1681
## published_after_2010 -0.2212 1.3848 -0.1597 0.8731 -2.9354 2.4930
##
## ---
## 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-yearsPrediction interval for reconviction
forensic_patient_reconviction$yi <- log((forensic_patient_reconviction$reconviction + 0.5) /
(forensic_patient_reconviction$total_year_complete + 0.5))
forensic_patient_reconviction$vi <- 1 / (forensic_patient_reconviction$reconviction + 0.5) +
1 / (forensic_patient_reconviction$total_year_complete + 0.5)
res_reconvict <- rma(yi = yi, vi = vi, data = forensic_patient_reconviction, method = "REML")
predict(res_reconvict, transf = exp) # Back-transforms to incidence rate scale##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.0314 0.0223 0.0442 0.0067 0.1465# Extract parameters
mu <- res_reconvict$b[1] # pooled log incidence rate
se <- res_reconvict$se[1] # standard error of pooled effect
tau2 <- res_reconvict$tau2 # between-study variance (tau squared)
# 95% CI of pooled effect (log scale)
ci_lb_log <- mu - 1.96 * se
ci_ub_log <- mu + 1.96 * se
# 95% Prediction interval (log scale)
pi_lb_log <- mu - 1.96 * sqrt(se^2 + tau2)
pi_ub_log <- mu + 1.96 * sqrt(se^2 + tau2)
# Back-transform to incidence rate scale
ci_lb <- exp(ci_lb_log)
ci_ub <- exp(ci_ub_log)
pi_lb <- exp(pi_lb_log)
pi_ub <- exp(pi_ub_log)
# Scale to per 100,000 person-years
ci_lb_100k <- ci_lb * 100000
ci_ub_100k <- ci_ub * 100000
pi_lb_100k <- pi_lb * 100000
pi_ub_100k <- pi_ub * 100000
pooled_rate_100k <- exp(mu) * 100000
# Print results
cat("Pooled incidence rate per 100,000 person-years:\n")## Pooled incidence rate per 100,000 person-years:cat(sprintf("%.2f [%.2f, %.2f]\n", pooled_rate_100k, ci_lb_100k, ci_ub_100k))## 3138.82 [2229.32, 4419.37]cat("\nPrediction interval per 100,000 person-years:\n")##
## Prediction interval per 100,000 person-years:cat(sprintf("[%.2f, %.2f]\n", pi_lb_100k, pi_ub_100k))## [672.41, 14652.08]## Borenstein's 95% prediction interval 672 - 14,650 per 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, 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)
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##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-yearsPrediction interval for readmission
forensic_patient_rehosp$yi <- log(forensic_patient_rehosp$readmission / forensic_patient_rehosp$total_year_complete)
forensic_patient_rehosp$vi <- 1 / forensic_patient_rehosp$readmission
res_rehosp <- rma(yi = yi, vi = vi, data = forensic_patient_rehosp, method = "REML")
predict(res_rehosp , transf = exp) # Back-transforms to incidence rate scale##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.0857 0.0602 0.1219 0.0163 0.4508# Extract parameters
mu <- res_rehosp$b[1] # pooled log incidence rate
se <- res_rehosp$se[1] # standard error of pooled effect
tau2 <- res_rehosp$tau2 # between-study variance (tau squared)
# 95% CI of pooled effect (log scale)
ci_lb_log <- mu - 1.96 * se
ci_ub_log <- mu + 1.96 * se
# 95% Prediction interval (log scale)
pi_lb_log <- mu - 1.96 * sqrt(se^2 + tau2)
pi_ub_log <- mu + 1.96 * sqrt(se^2 + tau2)
# Back-transform to incidence rate scale
ci_lb <- exp(ci_lb_log)
ci_ub <- exp(ci_ub_log)
pi_lb <- exp(pi_lb_log)
pi_ub <- exp(pi_ub_log)
# Scale to per 100,000 person-years
ci_lb_100k <- ci_lb * 100000
ci_ub_100k <- ci_ub * 100000
pi_lb_100k <- pi_lb * 100000
pi_ub_100k <- pi_ub * 100000
pooled_rate_100k <- exp(mu) * 100000
# Print results
cat("Pooled incidence rate per 100,000 person-years:\n")## Pooled incidence rate per 100,000 person-years:cat(sprintf("%.2f [%.2f, %.2f]\n", pooled_rate_100k, ci_lb_100k, ci_ub_100k))## 8567.57 [6022.66, 12187.87]cat("\nPrediction interval per 100,000 person-years:\n")##
## Prediction interval per 100,000 person-years:cat(sprintf("[%.2f, %.2f]\n", pi_lb_100k, pi_ub_100k))## [1628.33, 45078.90]## Borenstein's 95% prediction interval 1628 - 45078 per 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,
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)
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-yearsPrediction interval for all-cause mortality
forensic_patient_mortality$yi <- log(forensic_patient_mortality$all_cause_mortality / forensic_patient_mortality$total_year_complete)
forensic_patient_mortality$vi <- 1 / forensic_patient_mortality$all_cause_mortality
res_allcause <- rma(yi = yi, vi = vi, data = forensic_patient_mortality, method = "REML")
predict(res_allcause , transf = exp) # Back-transforms to incidence rate scale##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.0177 0.0127 0.0247 0.0069 0.0452# Extract parameters
mu <- res_allcause$b[1] # pooled log incidence rate
se <- res_allcause$se[1] # standard error of pooled effect
tau2 <- res_allcause$tau2 # between-study variance (tau squared)
# 95% CI of pooled effect (log scale)
ci_lb_log <- mu - 1.96 * se
ci_ub_log <- mu + 1.96 * se
# 95% Prediction interval (log scale)
pi_lb_log <- mu - 1.96 * sqrt(se^2 + tau2)
pi_ub_log <- mu + 1.96 * sqrt(se^2 + tau2)
# Back-transform to incidence rate scale
ci_lb <- exp(ci_lb_log)
ci_ub <- exp(ci_ub_log)
pi_lb <- exp(pi_lb_log)
pi_ub <- exp(pi_ub_log)
# Scale to per 100,000 person-years
ci_lb_100k <- ci_lb * 100000
ci_ub_100k <- ci_ub * 100000
pi_lb_100k <- pi_lb * 100000
pi_ub_100k <- pi_ub * 100000
pooled_rate_100k <- exp(mu) * 100000
# Print results
cat("Pooled incidence rate per 100,000 person-years:\n")## Pooled incidence rate per 100,000 person-years:cat(sprintf("%.2f [%.2f, %.2f]\n", pooled_rate_100k, ci_lb_100k, ci_ub_100k))## 1769.94 [1268.42, 2469.75]cat("\nPrediction interval per 100,000 person-years:\n")##
## Prediction interval per 100,000 person-years:cat(sprintf("[%.2f, %.2f]\n", pi_lb_100k, pi_ub_100k))## [693.64, 4516.30]## Borenstein's 95% prediction interval 693 - 4516 per 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,
xlim = c(-200, 5000))
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)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-yearsPrediction interval for suicide
forensic_patient_suicide$yi <- log(forensic_patient_suicide$suicide / forensic_patient_suicide$total_year_complete)
forensic_patient_suicide$vi <- 1 / forensic_patient_suicide$suicide
res_suicide <- rma(yi = yi, vi = vi, data = forensic_patient_suicide, method = "REML")
predict(res_suicide , transf = exp) # Back-transforms to incidence rate scale##
## pred ci.lb ci.ub pi.lb pi.ub
## 0.0039 0.0029 0.0053 0.0022 0.0069# Extract parameters
mu <- res_suicide$b[1] # pooled log incidence rate
se <- res_suicide$se[1] # standard error of pooled effect
tau2 <- res_suicide$tau2 # between-study variance (tau squared)
# 95% CI of pooled effect (log scale)
ci_lb_log <- mu - 1.96 * se
ci_ub_log <- mu + 1.96 * se
# 95% Prediction interval (log scale)
pi_lb_log <- mu - 1.96 * sqrt(se^2 + tau2)
pi_ub_log <- mu + 1.96 * sqrt(se^2 + tau2)
# Back-transform to incidence rate scale
ci_lb <- exp(ci_lb_log)
ci_ub <- exp(ci_ub_log)
pi_lb <- exp(pi_lb_log)
pi_ub <- exp(pi_ub_log)
# Scale to per 100,000 person-years
ci_lb_100k <- ci_lb * 100000
ci_ub_100k <- ci_ub * 100000
pi_lb_100k <- pi_lb * 100000
pi_ub_100k <- pi_ub * 100000
pooled_rate_100k <- exp(mu) * 100000
# Print results
cat("Pooled incidence rate per 100,000 person-years:\n")## Pooled incidence rate per 100,000 person-years:cat(sprintf("%.2f [%.2f, %.2f]\n", pooled_rate_100k, ci_lb_100k, ci_ub_100k))## 389.64 [288.98, 525.35]cat("\nPrediction interval per 100,000 person-years:\n")##
## Prediction interval per 100,000 person-years:cat(sprintf("[%.2f, %.2f]\n", pi_lb_100k, pi_ub_100k))## [219.82, 690.64]## Borenstein's 95% prediction interval 219 - 690 per 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,
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)