# 1. Bayesian meta-analysis with SMD as the outcome measure

# Set working directory and random seed, load packages, and define the model
library(dplyr)
setwd("D:/Rscript")
set.seed(12345)   # Set the seed
pacman::p_load(R2jags, tidyverse)

# Define the Bayesian model for SMD
model.smd=function(){
  for (i in 1:ns){
    prec[i] <- 1/v[i]
    y[i] ~ dnorm(delta[i], prec[i])
    delta[i] ~ dnorm(d, dprec)
  }
  d ~ dnorm(0, 1.0E-5)
  tau ~ dunif(0,2)
  tau.sq <- tau * tau
  dprec <- 1/(tau.sq)
}

# Create a dataset named PNFSMD.dat, after preprocessing to calculate SMD and its variance
PNFSMD.dat = tibble(
  nt = c(14,13,13,30,33),
  mt = c(-3.1,-2.51,-2.64,-5.4,-2.2),
  st = c(1.1,0.17,0.99,1.3,0.5),
  nc = c(13,14,16,30,33),
  mc = c(-0.1,-0.5,0.05,-3.5,-1.1),
  sc = c(1.1,0.95,0.08,1.7,0.9)
) %>%
  mutate(y = (1 - 3/(4*(nt+nc)-9)) * (mt - mc) / sqrt(((nt - 1)*st^2 + (nc - 1)*sc^2)/(nt + nc - 2))) %>%
  mutate(v = (nt + nc)/(nt*nc) + y^2/(2*(nt + nc - 3.94))) %>%
  select(y, v)

# Load the data, set initial values and parameters to monitor, fit the model using jags(), and display results
data = list(y = PNFSMD.dat$y, v = PNFSMD.dat$v, ns = length(PNFSMD.dat$y))
inits = function(){ list(d = 0) }
paras = c("d", "tau", "tau.sq")

res.smd = jags(data, inits, paras, model.file = model.smd, n.iter = 50000)
print(res.smd)

#  2. Kappa 
#2.1.  Simple kappa consistency test for unordered variables
data1 <- matrix(
  c(9,1,1,913),
  2,2,
  dimnames = list(
    R2 = c("Included","Excluded"),
    R1 = c("Included","Excluded")
  )
)
data1

install.packages("fmsb")
library(fmsb)

easykappa <- Kappa.test(data1)
easykappa

#2.2. Weighted kappa consistency test for ordered variables
# The matrix is entered column by column from left to right
data2 <- matrix(
  c(11,0,0, 1,13,1, 0,0,4),
  3,3,
  dimnames = list(
    R2 = c("Low","Medium","High"),
    R1 = c("Low","Medium","High")
  )
)
data2

# Function to convert frequency tables to case-level data
counts_to_cases <- function(x, count.col = "Freq") {
  if(!inherits(x, "table")) x <- as.table(as.matrix(x))
  x <- as.data.frame(x)
  # Generate row indices to replicate rows
  idx <- rep.int(seq_len(nrow(x)), x[[count.col]])
  # Remove the frequency column
  x[[count.col]] <- NULL
  # Create the expanded dataset
  x <- x[idx, ]
  rownames(x) <- 1:nrow(x)
  x
}

# Convert frequency table to individual cases
data2a <- counts_to_cases(x = data2, count.col = "Freq")
data2a

# Convert factors to numeric
mat_factor2 <- factor(data2a$R2, levels = c("Low", "Medium", "High"))
mat_factor1 <- factor(data2a$R1, levels = c("Low", "Medium", "High"))
mat_num1 <- as.numeric(mat_factor1)
mat_num2 <- as.numeric(mat_factor2)
datakappa <- cbind(mat_num1, mat_num2)

# Calculate the weighted kappa
install.packages("psy")
library(psy)
wkappa(datakappa, "absolute")

# Calculate the 95% confidence interval for kappa
install.packages("boot")
library(boot)

wkappa.boot <- function(data, x) {
  wkappa(data[x,])[[3]]
}

res <- boot(datakappa, wkappa.boot, 1000)
quantile(res$t, c(0.025, 0.975))
boot.ci(res, type = "bca")