######################################
##### Simple Regression Commands  ####
#####       Paris 8 2019          ####
######################################

## Packages
## (MASS) -  comes pre-installed
## (rms) -   download / install with package manager

## There are two packages for doing a logistic regression analysis.
## We use both glm() function (package - MASS) + the lrm() function (package rms). 
## The glm() function is simpler, but gives you less information. We begin with glm.

#################
##  Load Data  ##
#################

# Data format - make sure it is a flat data frame, not a numerical summary / cross-tabulations of the results.
# Load the entire data frame, we will choose which variables to use "in" R.

df <- read.table(file.choose(), header=T, sep="\t")

####################
##  Examine Data  ##
####################

# To look at the data and check for spelling mistakes and to make sure all has worked

summary(df)

# Some of your variables will have more then 6 levels (features / categories)
# R only lists the 7 most common features. To see the rest, to check for spelling issues and small counts, use the function table()
# In following line, replace "Name_of_Variable" with the name of the variable (factor / column name) with more than 6 levels
# Copy and paste from the output of the summary() command works well

table(df $ Name_of_Variable)

# Make sure you clean the data (spelling, small counts etc.) in the original data - yr spreadsheet (e.g. Excel) or database (e.g. Filemaker)
# When the data are 'clean', reload with the read.table() function above.

#####################################
#####  GLM Logistic Regression  #####
#####################################

# Use the output from the summary() function to choose (and copy and paste) the response variable and the predictor variables
# Add all predictor variables as you believe may have a role in the outcome, separating each by a '+' sign.
# If you get an error message saying the model cannot converge: there are two possible reasons
# -  you have a variable that has features with only a few occurrences (12 is minimum, but aim for 20)
# -  you have two variables that have a feature which co-occurs perfectly(predicts in the same way). Such variables cannot be submitted simultaneously. See multicollinearity below.

# Windows users: You cannot select the 'predict_variable' in the command line below, you must use back space and delete the 'predictor_variable' before you add the 'real' variable names
# Mac users: simply 'copy' the predictor variable name, then highlight 'predict_variable' in the line below and 'paste'

M = glm(DV ~ IV + IV + IV, data = df, family = "binomial")

summary(M)

#####################################
#####   LRM Logistic Regression  ####
#####################################

# lrm does the same as glm, but it does not offer stars to show significance and you cannot plot the results. 
# The table of co-efficient should be the same as a glm (save it rounds down very small numbers)
# lrm produces the model statistics, which allow you to determine the predictive strength of the model
# Don't forget you need to load the libraries - lrm() requires rms (see above)

library(rms)

LRM = lrm(Response_Variable ~ predict_variable + predict_variable + predict_variable, data = df, x=T, y=T)
LRM

#########################
##  Model Diagnostics  ##
#########################

# Once you have 1 or 2 models that are parsimonious and interpretable, we need to check the model for a range of problems
# The most common problems are: multicollinearity, over-dispersion, over-fitting
# The following descriptions and R-commands will help you check and improve your models
# These steps take much more time than actually fitting (doing) a logistic regression
# Model diagnostics are not restricted to what is below, but this is a good place to start

## Multicollinearity ###

# This concept is crucial to logistic regression and is typically the biggest problem in corpus linguistics 
# Logistic regression does NOT assume that the data is normally distributed etc, but it does assume that a model is 'orthogonal'
# Each predictor variable must predict the outcome in a different way. If you have two variables that are 'similar', then this artificially improves the results (possible Type 1 error)
# Having no multicollinearity / having an orthogonal model is essential

vif(ModGLM)

# The vif() function calculates the 'variance inflation factor'
# This is a figure that tells you if a given feature (level) is collinear in the model. 
# There is a debate over what is an acceptable level of multicollinearity. 
# Any number at or above 10 is clearly a problem, and the to which it belongs must removed.
# Any number above 4 should be reported, as many (if not most) people consider 4 to be the upper-limit of acceptability.
# Some pople consider 2.5 to be the maximum VIF.
# This should work on both glm and lrm models

## Bootstrap Validation ##

# Bootstrapping is a randomisation method which allows you estimate probabilities and significance.
# We can apply it to the model to 'check' the model statistics, such as the Pseudo R2. 
# The following is taken from Baayen (2008: 205f).

validate(ModLRM, bw = T, B = 200)

# If the bootstrapped R2 is considerably lower than that produced with the lrm() function, then the model suffers from over-dispersion and its predictive strength should be interpreted with great care.

#########################
## The Basic Commands ###
#########################

#################################
###  Logistic Regression GLM ####
#################################

# Logistic Regression needs its data in data frame format

M <- glm(Response_Variable ~ Variable1 + Variable2 + Variable3, data = df, family = "binomial")
summary(M)

###################################
###  Logistic Regression in LRM ###
###################################

# Logistic Regression needs its data in data frame format

LRM <- lrm(Response_Variable ~  Variable1 + Variable2 + Variable3, data = df, x=T, y=T)
LRM