# Edsp 589
# Fall 2021
# c.j.anderson
#
# Function for range-of-influence statistics for detection of 
#  possibly miss-classified observations on the response/outcome
#  variable...with ESR examples.
#
# The function find ROI for regression coefficients and deviance.
#  This can be changed without too much difficulty
#

########################################################################
# roi function(indata, model.formula)                                  #                   
########################################################################
# where indata= your unalterated data set-up
#       model.formula you want to use (see examples below on how to create a formula object).
#
# Output:
#  [[1]]  The model with unaltered data
#  [[2]]  The coefficents from each model
#  [[3]]  ROI statistics for coefficents and last column is for deviance
#
#
# set y < your.response.variable
#
#
#########################################

roi <- function(indata,model.formula) {                                # SET-UP
        copy.data <- indata                                            # make a copy of data to use
        n <- length(copy.data$y)                                       # number of observations
	copy.data$id <- seq(1:n)                                       # add an id variable 
        logit0 <- glm(model.formula,data=copy.data,family=binomial)    # model with no changes to data
        original.coeff <- logit0$coefficients                          # coefficients from this model
        original.dev <- logit0$deviance 		               # deviance from this model
	n.coeff <- length(logit0$coefficients)                 	       # number of coefficients
	coeff.saved <- matrix(99,nrow=n,ncol=(n.coeff+1) )             # the +1 is for deviance
        roi <- matrix(99,nrow=n,ncol=(n.coeff+1)   )                   # roi will be saved here                                       
 
        for (loop in 1:n) {                                            # LOOP THROUGH CASES
 		   workdata <- copy.data                                # workdata will be one with changed y
	       if (workdata$id[loop]==loop) {                           # change 0 to 1 or 1 to 0
		      workdata$y[loop] = 1 - workdata$y[loop]
		   }
		   logit1 <- glm(model.formula,data=workdata,family=binomial)  # fit model to altered data
		   coeff.saved[loop,] <- c(logit1$coefficients,logit1$deviance)# save coefficients and deviance
		   roi[loop,1:n.coeff] <- logit1$coefficients - original.coeff # roi for coefficients
		   roi[loop,n.coeff+1] <- logit1$deviance - original.dev       # roi for deviance		   
		}                           
		outlist <- list(logit0,coeff.saved,roi)                        # objects output       
        return(outlist)
}		
###########################################################################
# End of function                                                         #
###########################################################################

###########################################################################
# Examples:  Use esr data
###########################################################################
setwd("C:\\Users\\cja\\Dropbox\\edps 589\\5 logistic regression")
esr <- read.table("esr_data.txt",header=TRUE)
#
#
# Need to make your response/outcome variable
esr$y <- esr$response                # Need to make your response/outcome variable y

#
# Example 1
#

model.to.fit <- formula(y ~ fibrin)  # this is the model you want to check

roi.list <- roi(esr,model.to.fit)    # output is a matrix object 


# Post run looking
roi.stat <- as.data.frame(roi.list[[3]])  # Make roi statistics a data frame
roi.stat$id <- seq(1:length(esr$response))# Add id to roi data frame so can id
                                          # specific cases.
plot(roi.stat$V2,roi.stat$V1, type="n",   # Pot data with points being id
     main="ROI: slope x intercept",       # and reference lines at 0,0
     xlab="Slope for Fibrin",
     ylab="Intercept", 
     col="blue")
text(roi.stat$V2,roi.stat$V1, roi.stat$id)
lines(c(0,0),c(min(roi.stat$V1),max(roi.stat$V1)),lty=2 )
lines(c(min(roi.stat$V2),max(roi.stat$V1)),c(0,0),lty=2 )


#
# Example 2
#
model2.to.fit <- formula(y ~ fibrin + globulin)  # this is the model you want to check

roi.list <- roi(esr,model2.to.fit)

roi.stat <- as.data.frame(roi.list[[3]])
roi.stat <- roi.stat[order(roi.stat$V4),]  # order by deviance
roi.stat