## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7, 
  fig.height = 4
)

## ----setup--------------------------------------------------------------------
library(CLRtools)
library(survival)
library(dplyr)

## -----------------------------------------------------------------------------
# Predictor variables to evaluate
var_to_test <- c('height','weight', 'bmi', 'priorfrac', 'premeno', 'momfrac', 'armassist','smoke', 'raterisk')

# Define OR scaling factors (used to interpret ORs per meaningful unit)
OR_values <- c(10, 10, 3, 1, 1, 1, 1, 1, 1, 1)

# Generate the univariable models table
univariable.clogmodels(glow11m, yval = 'fracture',xval = var_to_test, strata='pair', OR=TRUE, inc.or = OR_values)

# Obtaining the discordant pairs
discordant.pairs(glow11m, outcome = 'fracture', strata = 'pair')

## -----------------------------------------------------------------------------
summary(clogit(as.numeric(fracture) ~ height + weight + bmi + priorfrac + premeno + momfrac + armassist + raterisk + strata(pair), data = glow11m))

## -----------------------------------------------------------------------------
# Variables that were not significant in Step 2 to check
candidate.exclusion <- c('premeno', 'raterisk') 

# Significant variables in Step 2
var.preliminar <- c('height', 'weight', 'bmi','priorfrac',  'momfrac', 'armassist') 

# Check for confounding 
check_coef_change(data = glow11m, yval = 'fracture', xpre = var.preliminar, xcheck = candidate.exclusion, strata = 'pair')

## -----------------------------------------------------------------------------
summary(clogit(as.numeric(fracture) ~ weight + bmi + priorfrac + momfrac + armassist + strata(pair), data = glow11m))

## -----------------------------------------------------------------------------
# Variables excluded in Step 1 to check
excluded<-c('smoke')

# Preliminary model variables (retained after Step 2 and 3)
var.preliminar <- c('weight','bmi', 'priorfrac', 'momfrac','armassist') 

check_coef_significant(data = glow11m, yval = 'fracture', xpre = var.preliminar, xcheck = excluded, strata = 'pair')

## -----------------------------------------------------------------------------
var.maineffects<-c('weight', 'bmi','priorfrac','momfrac','armassist')

check_interactions(data = glow11m, variables = var.maineffects, yval = 'fracture', rounding = 4, strata = 'pair')

## -----------------------------------------------------------------------------
final.model <- clogit(as.numeric(fracture) ~ weight + bmi + priorfrac + momfrac + armassist + strata(pair), data = glow11m)
summary(final.model)

## -----------------------------------------------------------------------------
# Converting the outcome from "Yes"/"No" to binary 0/1, as required by the function
glow11m$fracture <- ifelse(glow11m$fracture == "Yes", 1, 0)

# 
head(residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$data.results)
residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$leverage
residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$change.Pearsonchi
residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$cooks
residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$m.Pearsonchi
residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$mcooks

## -----------------------------------------------------------------------------
# Obtain residual diagnostics for the conditional logistic regression model.
residuals.data<-residuals_clog(final.model, strata = glow11m$pair, y = glow11m$fracture, plot = TRUE)$data.results

# Summarize the delta chi-squared and delta beta for each matched pair.
glow.match<-residuals.data %>% group_by(pair) %>% dplyr::summarise(sum_deltachi=sum(deltachi), sum_deltabeta=sum(deltabeta))

# Identify individual observations considered potential outliers based on thresholds from the plot
outliers<-subset(residuals.data, deltachi>3 | deltabeta>0.15 )
outliers.match<-subset(glow.match, sum_deltachi>3 | sum_deltabeta >0.15)

# Subset all observations belonging to the identified outlier pairs.
glow.outieliers<-subset(residuals.data, pair %in% outliers$pair)

# Reorder and display the final table of potential outliers and their diagnostics.
glow.outieliers.t<-glow.outieliers[,c(7,1:5,8,12,13,10)]
glow.outieliers.t
outliers.match