## Replication script for examples used in RESI: An R Package for Robust Effect Sizes
## Megan Jones, Kaidi Kang, Simon Vandekar

# install.packages("RESI")
# install.packages("sars")
# install_version("RESI")
# install.packages("survival")

# RESI on linear model
library(RESI)
mod_lm <- lm(charges ~ region * age + sex + bmi, data = insurance)
summary(mod_lm)

set.seed(0827)
resi_obj_lm <- resi(mod_lm)
summary(resi_obj_lm)

anova(resi_obj_lm)

omnibus(resi_obj_lm)

library(ggplot2)
plot(resi_obj_lm)
ggplot(anova(resi_obj_lm))

set.seed(0827)
resi_obj_lm2 <- resi(mod_lm, vcov.args = list(type = "HC0"), Anova.args = list(type = 3))
resi_obj_lm2

# RESI on nonlinear model
data(niering, package = "sars")
head(niering)

# run to display an error
# mod_nls <- nls(s ~ c*a^z, data = niering, start = list(c = 2, z = 0.5))

mod_nls <- nls(s ~ c*a^z, data = niering, start = list(c = 3, z = 0.25))
summary(mod_nls)

set.seed(0827)
resi_obj_nls <- resi(mod_nls, data = niering, boot.method = "bayes")
resi_obj_nls

mod_nls2 <- nls(s ~ c*a^z, data = niering, start = list(c = coef(mod_nls)[1],
                                                        z = coef(mod_nls)[2]))
set.seed(0827)
resi(mod_nls2, data = niering, boot.method = "bayes")

# RESI on survival model
library(survival)
set.seed(0827)
mod_surv <- survreg(Surv(time, status) ~ age + sex + ph.karno,
                    data = survival::lung, dist="weibull", robust = TRUE)
mod_surv_reduced <- survreg(Surv(time, status) ~ ph.karno,
                            data = survival::lung, dist="weibull", robust = TRUE)
resi_obj_surv <- resi(mod_surv, mod_surv_reduced, data = survival::lung,
                      unbiased = FALSE, store.boot = TRUE, alpha = c(0.05, 0.1),
                      nboot = 1500)
resi_obj_surv

set.seed(0827)
omnibus(resi(mod_surv, data = survival::lung,
     unbiased = FALSE, alpha = c(0.05, 0.1), nboot = 1500))

summary(resi_obj_surv, alpha = c(0.001, 0.01))

anova(resi_obj_surv, alpha = c(0.001, 0.01))


sessionInfo()