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

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

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

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

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

## -----------------------------------------------------------------------------
summary(glm(fracture ~ age + height + priorfrac + momfrac + armassist + raterisk, family = binomial, data = glow500))

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

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

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

## -----------------------------------------------------------------------------
# Variables excluded in Step 1 to check
excluded<-c('weight', 'bmi', 'premeno','smoke')

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

check_coef_significant(data = glow500, yval = 'fracture', xpre = var.preliminar, xcheck = excluded)

## -----------------------------------------------------------------------------
# Transforming raterisk
glow500<-glow500 %>% mutate(raterisk_cat=case_when((raterisk=='Less'| raterisk=='Same')~'C1',
                                                   raterisk=='Greater'~'C2'))

summary(glm(fracture ~ age + height + priorfrac + momfrac + armassist + raterisk_cat, family = binomial, data = glow500))

## -----------------------------------------------------------------------------
var.maineffects<-c('age', 'height', 'priorfrac', 'momfrac', 'armassist', 'raterisk_cat')

check_interactions(data = glow500, variables = var.maineffects, yval = 'fracture', rounding = 4)

## -----------------------------------------------------------------------------
final.model <- glm(fracture ~ age + height + priorfrac + momfrac + armassist + raterisk_cat + age*priorfrac + momfrac*armassist, family = binomial, data = glow500)
summary(final.model)

## -----------------------------------------------------------------------------
DRtest(final.model, g = 10)

## -----------------------------------------------------------------------------
osius_rojek(final.model)

## -----------------------------------------------------------------------------
stukels_test(final.model) 

## -----------------------------------------------------------------------------
cutpoints(final.model, cmin = 0.05, cmax = 0.75, byval = 0.05, plot = TRUE)

## -----------------------------------------------------------------------------
diagnosticplots_class(final.model)

## -----------------------------------------------------------------------------
# Measures assuming J = n, with J being the number of covariate patterns
r_measures(final.model)

## -----------------------------------------------------------------------------
# Measures assuming J < n
rcv_measures(final.model)

## -----------------------------------------------------------------------------
head(residuals_logistic(final.model)$x.cv)
residuals_logistic(final.model)$leverage
residuals_logistic(final.model)$change.Pearsonchi
residuals_logistic(final.model)$change.deviance
residuals_logistic(final.model)$cooks
residuals_logistic(final.model)$change.Pb

## -----------------------------------------------------------------------------
# Extract the residuals and other measures
residuals.data<-residuals_logistic(final.model)$x.cv

# Restore variable names for readability
colnames(residuals.data)[c(4:7)]<-c('priorfrac','momfrac','armassist','raterisk_cat')

# Flag observations with large influence or poor fit, and filter them
residuals.data<-residuals.data%>% mutate(outliers=ifelse(delta.chi>10 | delta.deviance >4.5 | delta.beta>0.12, 1,0))
out<-residuals.data %>% filter(outliers==1)

# Select and format relevant columns for display
out<-out[,c(2:7, 10:12, 18:19, 13, 17)]
out[order(out[, 1], decreasing = FALSE), ]