# the library() function loads the modules that we need
library(gmodels)
# next, we open the data file we want to analyze and add it to the search string
Admit <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
attach(Admit)
# in this data, we'll use three predictors (GRE, GPA, and rank) to explain variation
# in admission status (admit, not admit)
# first, let's look at the variables
summary(Admit)
# let's look at the categorical predictor and the criterion in a crosstable to make
# sure that there aren't empty cells.
CrossTable(admit, rank, prop.r=F, prop.c=T, prop.t=F, prop.chisq=F, chisq = F)
# since there are no empty cells, we can use this variable as a predictor
# but first we have to make a series of dummy vectors (binary variables)
# to identify each level of rank. I use the factor() command and create a
# new variable, rankf. When we use rankf in the regression model, it will use
# the dummy vectors automatically.
rankf <- factor(rank)
# now our logistic regression model
summary(glm(admit ~ gre + gpa + rankf, family=binomial(), data=Admit))
# let's take a look at another model using out USCounties data
USCounties <- read.csv("http://www.courseserve.info/files/USCounties.csv")
attach(USCounties)
# first we need to fix the income per capita data
CA05N0030_05[CA05N0030_05==-3000000000]<-NA
summary(glm(LowEduc04 ~ PopGr00_05 + IntlMig05 + MedianAge05 + PctWhite05 + CA05N0030_05, family=binomial(), data=USCounties))
exp(coef(glm(LowEduc04 ~ PopGr00_05 + IntlMig05 + MedianAge05 + PctWhite05 + CA05N0030_05, family=binomial(), data=USCounties)))
logLik(glm(LowEduc04 ~ PopGr00_05 + IntlMig05 + MedianAge05 + PctWhite05 + CA05N0030_05, family=binomial(), data=USCounties))
# to test the model for significance, we use the anova() function with Chi-squared
# we're comparing the added explanatory power of the predictors compared to the empty model (only the intercept)
anova(glm(LowEduc04 ~ 1, family=binomial(), data=USCounties), glm(LowEduc04 ~ PopGr00_05 + IntlMig05 + MedianAge05 + PctWhite05 + CA05N0030_05, family=binomial(), data=USCounties), test="Chisq")