---
title: 'Multiple Logstic Regression: HSB'
author: "Carolyn J. Anderson"
date: '2022-10-17'
output: 
   html_document:
     toc: yes
---

These are analyses reported in lecture notes on logistic regression.  

# Set up

## libraries

```{r libs}}
library(vcd)
library(ResourceSelection)
library(pROC)               
library(DescTools)
library(texreg)
```

## Data

### In long format (data frame):
```{r long_df}
setwd("C:/Users/cja/Dropbox/edps 589/5 logistic regression")
hsb <- read.csv(file="hsb_data.csv", header=TRUE, sep=",")
names(hsb)
head(hsb)
```

### Long table:
```{r long_table}
var.levels <- expand.grid(academic=c("no", "yes"),
                          sctyp=c("public","private"),
			  ses = c("low","middle","high"))
( hsb3way <- data.frame(var.levels, count=c(91, 40,  4,  4, 138, 111,
                                            14, 36, 44, 82,   1,  35 )))				  
```

### Wide table:

```{r wide_table}
var.levels2 <- expand.grid(sctyp=c("public","private"),
			  ses = c("low","middle","high"))
hsb3way.no <- data.frame(var.levels2, no=c(40, 4, 111,  36, 82,   35 ))	 					  
hsb3way.yes <- data.frame(var.levels2, yes=c(91,4,136,14,44,1))					  
hsb3way.wide <- merge(hsb3way.yes,hsb3way.no, by=c("sctyp","ses"))
hsb3way.wide$ncases <- hsb3way.wide$yes + hsb3way.wide$no
(hsb3way.wide)
```

# Some Logit Models

## Independence
```{r logit_independence}
summary(  i.mod <- glm(cbind(yes,no) ~ 1 , data=hsb3way.wide, family=binomial))
```

## Homogeneous Associaiton

or all 2-way interactions
```{r logit_2way}
summary(  h.logit <- glm(cbind(yes,no)  ~ sctyp + ses , data=hsb3way.wide,family=binomial)  )

(X2 <- sum(residuals(h.logit,type="pearson")^2) )
```

# Logit Models as poisson regression models

```{r logit_as_loglinear}
# As a poisson --- one way to get  X2 goodness-of-fit statistics
summary( h.poi <- glm(count ~ sctyp + academic + ses + sctyp*academic + sctyp*ses + ses*academic, data=hsb3way,family=poisson) )
(X2 <- sum(residuals(h.poi,type="pearson")^2))

# This will give both G2 and X2
library(MASS)
summary( h.poi2 <- loglm(count ~ sctyp + academic + ses + sctyp*academic + sctyp*ses + ses*academic, data=hsb3way) )
```

# Identification Constraints

Different programs may use different ID constraints.  Different parameterizations on impact the value of the regression coefficients but not the fitted values, odds ration estimates, residuals, goodness-of-model fit to data, etc.

## Dummy: first set to zer0
The first is set to 0,
```{r default}
hsb3way.wide$low   <- ifelse(hsb3way.wide$ses=="low",1,0)
hsb3way.wide$mid   <- ifelse(hsb3way.wide$ses=="middle",1,0)
hsb3way.wide$hig   <- ifelse(hsb3way.wide$ses=="high",1,0)

summary( logit.bin <- glm(cbind(yes,no) ~  mid + hig + sctyp, family = binomial, data = hsb3way.wide) )

```

## Alternate Dummy Codes
By defaults glm in R will set the 1st coefficient to 0, but SAS sets the last to 0.  Below we switch to last coefficients

```{r alternate_dummy}
# -- change how Dummy codes are defined
hsb3way.wide$Dsctyp <- ifelse(hsb3way.wide$sctyp=="public",1,0)
hsb3way.wide$Dlow   <- ifelse(hsb3way.wide$ses=="low",1,0)
hsb3way.wide$Dmid   <- ifelse(hsb3way.wide$ses=="middle",1,0)

summary( logit.bin2 <- glm(cbind(yes,no) ~  Dlow + Dmid + Dsctyp, family = binomial, data = hsb3way.wide) )
```

## Effects codes

```{r effect_codes}
hsb3way.wide$Esctyp <- ifelse(hsb3way.wide$sctyp=="public",1,-1)
hsb3way.wide$Elow   <- ifelse(hsb3way.wide$ses=="low",1,0)
hsb3way.wide$Elow   <- ifelse(hsb3way.wide$ses=="high",-1, hsb3way.wide$Elow)
hsb3way.wide$Emid   <- ifelse(hsb3way.wide$ses=="middle",1,0)
hsb3way.wide$Emid   <- ifelse(hsb3way.wide$ses=="high",-1,hsb3way.wide$Emid)

summary( logit.bin3 <- glm(cbind(yes,no) ~  Elow + Emid + Esctyp, family = binomial, data = hsb3way.wide) )
```
## Comparison

```{r}
screenreg(list(logit.bin, logit.bin2, logit.bin3),
          custom.model.names=c("First=0", "Last =0", "Effect"))
```
Same model, just different value of the coefficients.

## Same odds ratios

Holding sctype constant:  public vs private odds  ratios

```{r odds_ratios}
(odds.dummy1 <- exp(0 - logit.bin$coefficients[4]))

(odds.dummy2 <-  exp(logit.bin2$coefficients[4]- 0))
 
(odds.effect3 <-  exp(logit.bin3$coefficients[4]+logit.bin3$coefficients[4]))
```

You can re-code school type to get exactly what are in the lecture notes

# Nominal versus Ordinal Explanatory Variable?

SES is definitely ordinal, but first treat is as a factor 
```{r ses_factor}
hsb3way.wide$xses  <- as.numeric(hsb3way.wide$ses)

# -- ses as factor
summary( logit.factor <- glm(cbind(yes,no) ~  ses + sctyp, family = binomial, data = hsb3way.wide) )
```

Note SES will be numeric (this is special case of when treated as a factor)
```{r ses_numeric}
summary( logit.numeric <- glm(cbind(yes,no) ~  xses + sctyp, family = binomial, data = hsb3way.wide) )
```

Compute goodness-of-fit statistics for logit.numeric
```{r}
hsb3way.wide$pi.numeric <- logit.numeric$fitted
Xsq.yes <- (hsb3way.wide$yes -     hsb3way.wide$pi.numeric*hsb3way.wide$ncases)**2   /(hsb3way.wide$pi.numeric*hsb3way.wide$ncases)
Xsq.no <-  (hsb3way.wide$no   - (1-hsb3way.wide$pi.numeric)*hsb3way.wide$ncases)**2 /((1-hsb3way.wide$pi.numeric)*hsb3way.wide$ncases)
(Xsq <- sum(Xsq.yes + Xsq.no))

Gsq.yes <- 2*hsb3way.wide$yes *log(hsb3way.wide$yes/(hsb3way.wide$ncases*hsb3way.wide$pi.numeric))
Gsq.no <- 2*hsb3way.wide$no * log(hsb3way.wide$no/(hsb3way.wide$ncases*(1-hsb3way.wide$pi.numeric)))
( Gsq <- sum(Gsq.yes+Gsq.no) )

# -- LR test of restriction on values for SES
anova(logit.numeric, logit.factor)
```

Can treat SES as numeric with evenly spaced scores.

## Odds Ratio with ordinal/numeric SES

Just some here
```{r some_odds_ordinal}
logit.numeric$coefficients

# -- Holding schtyp constant, odds of academic are 
exp(logit.numeric$coefficients[2])
# the odds for an increase of SES by 1 unit

# -- odds academic of low vs high ses
exp((3-1)*logit.numeric$coefficients[2])
```


# Fitting Data at Individual Level

## (1) Simple:  program ~ math
```{r model_1}
summary(logit.mod1 <- glm(program ~ math, data=hsb,family = binomial)) 
j <- order(hsb$math)	                                
plot(hsb$math[j],logit.mod1$fitted[j], 
    type='l',
	lwd=2,
	ylim=c(0,1),
	col="blue",
	main="Model 1: program ~ math")
```

## (2) SES as nominal
```{r model_2}
hsb$ses <- as.factor(hsb$ses)
summary(logit.mod2 <- glm(program ~ math + ses, data=hsb,family = binomial)) 

# - subset for graphing
hsb$fitted2 <- logit.mod2$fitted
low <- hsb[ which(hsb$ses=="1"),]
mid <- hsb[ which(hsb$ses=="2"),]
hi  <- hsb[ which(hsb$ses=="3"),]
j <- order(low$math)	                                
plot(low$math[j],low$fitted2[j], 
    type='l',
	lwd=2,
	ylim=c(0,1),
	xlim=c(30,80),
	col="blue",
	main="Model 2: program ~ math + ses (as factor) "
	)
j <- order(mid$math)	                                
lines(mid$math[j],mid$fitted2[j],col="green",lwd=2)
j <- order(hi$math)	                                
lines(hi$math[j],hi$fitted2[j],col="red",lwd=2)
legend("topleft",legend=c("high","middle","low"), 
       col=c("red","green","blue"), lty=1,lwd=2, cex=1)
```

## (3) numeric SES
```{r model_3}
hsb$ses <- as.numeric(hsb$ses)
summary(logit.mod3 <- glm(program ~ math + ses, data=hsb,family = binomial)) 

# - subset for graphing
hsb$fitted3 <- logit.mod3$fitted
low <- hsb[ which(hsb$ses=="1"),]
mid <- hsb[ which(hsb$ses=="2"),]
hi  <- hsb[ which(hsb$ses=="3"),]
j <- order(low$math)	                                
plot(low$math[j],low$fitted3[j], 
    type='l',
	lwd=2,
	ylim=c(0,1),
	xlim=c(30,80),
	col="blue",
	main="Model 3: program ~ math + ses (numeric) "
	)
j <- order(mid$math)	                                
lines(mid$math[j],mid$fitted3[j],col="green",lwd=2)
j <- order(hi$math)	                                
lines(hi$math[j],hi$fitted3[j],col="red",lwd=2)
legend("topleft",legend=c("high","middle","low"), 
       col=c("red","green","blue"), lty=1,lwd=2, cex=1)
```

## (4) program ~ math + ses + sctyp
```{r model_4}
hsb$sctyp <- as.factor(hsb$sctyp)
summary(logit.mod4 <- glm(program ~ math + ses + sctyp, data=hsb,family = binomial)) 

# - subset for graphing
hsb$fitted4 <- logit.mod4$fitted
low.public  <- hsb[ (which(hsb$ses==1 & hsb$sctyp=="1")),]
low.private <- hsb[ (which(hsb$ses==1 & hsb$sctyp=="2")),]
mid.public  <- hsb[ (which(hsb$ses==2 & hsb$sctyp=="1")),]
mid.private <- hsb[ (which(hsb$ses==2 & hsb$sctyp=="2")),]
hi.public   <- hsb[ (which(hsb$ses==3 & hsb$sctyp=="1")),]
hi.private  <- hsb[ (which(hsb$ses==3 & hsb$sctyp=="2")),]
j <- order(low.public$math)	                               
plot(low.public$math[j],low.public$fitted4[j], 
    type='l',
	lwd=2,
	ylim=c(0,1),
	xlim=c(30,80),
	col="blue",
	main="Model 4: program ~ math + ses (numeric) + school type  "
	)
j <- order(low.private$math)	                                
lines(low.private$math[j],low.private$fitted4[j],col="darkblue",lwd=2)
j <- order(mid.public$math)	                                
lines(mid.public$math[j],mid.public$fitted4[j],col="green",lwd=2)
j <- order(mid.private$math)	                                
lines(mid.private$math[j],mid.private$fitted4[j],col="darkgreen",lwd=2)
j <- order(hi.public$math)	                                
lines(hi.public$math[j],hi.public$fitted4[j],col="red",lwd=2)
j <- order(hi.private$math)	                                
lines(hi.private$math[j],hi.private$fitted4[j],col="darkred",lwd=2)
legend("bottomright",legend=c("private high SES", "private middle SES","private low SES","public high SES", "public middle SES", "public low SES "), 
       col=c("darkred","darkgreen","darkblue","red","green","blue"), lty=1,lwd=2, cex=1)
```

## (5) An interaction  program ~ math + ses + sctyp + ses*sctyp
```{r model_5}
hsb$sctyp <- as.factor(hsb$sctyp)
summary(logit.mod5 <- glm(program ~ math + ses + sctyp + ses*sctyp, data=hsb,family = binomial)) 

# - subset for graphing
hsb$fitted5 <- logit.mod5$fitted
low.public  <- hsb[ (which(hsb$ses==1 & hsb$sctyp=="1")),]
low.private <- hsb[ (which(hsb$ses==1 & hsb$sctyp=="2")),]
mid.public  <- hsb[ (which(hsb$ses==2 & hsb$sctyp=="1")),]
mid.private <- hsb[ (which(hsb$ses==2 & hsb$sctyp=="2")),]
hi.public   <- hsb[ (which(hsb$ses==3 & hsb$sctyp=="1")),]
hi.private  <- hsb[ (which(hsb$ses==3 & hsb$sctyp=="2")),]
j <- order(low.public$math)	                               
plot(low.public$math[j],low.public$fitted5[j], 
    type='l',
	lwd=2,
	ylim=c(0,1),
	xlim=c(30,80),
	col="blue",
	main="Model 5: program ~ math + ses + schyp + ses*schtyp "
	)
j <- order(low.private$math)	                                
lines(low.private$math[j],low.private$fitted5[j],col="darkblue",lwd=2)
j <- order(mid.public$math)	                                
lines(mid.public$math[j],mid.public$fitted5[j],col="green",lwd=2)
j <- order(mid.private$math)	                                
lines(mid.private$math[j],mid.private$fitted5[j],col="darkgreen",lwd=2)
j <- order(hi.public$math)	                                
lines(hi.public$math[j],hi.public$fitted5[j],col="red",lwd=2)
j <- order(hi.private$math)	                                
lines(hi.private$math[j],hi.private$fitted5[j],col="darkred",lwd=2)

legend("bottomright",legend=c("private high SES", "private middle SES","private low SES","public high SES", "public middle SES", "public low SES "), 
       col=c("darkred","darkgreen","darkblue","red","green","blue"), lty=1,lwd=2, cex=1)
```

## Diagonistics for model 5

### Hosmer-Lemeshow 
```{r hosmwer_lemshow}
hoslem.test(logit.mod5$y, fitted(logit.mod5), g=10)
```


### qq plot residuals

Using Pearson (adjusted standardized residuals
```{r qqadj}
resid <- rstandard(logit.mod5, pearson=TRUE)
qqnorm(resid,
       main="QQ plot of  Pearson  Residuals",
	     col="blue")
```
### ROC and concordance
#### Concordance index
```{r auc}
Cstat(logit.mod1,resp=hsb$program)
Cstat(logit.mod2,resp=hsb$program)
Cstat(logit.mod3,resp=hsb$program)
Cstat(logit.mod4,resp=hsb$program)
Cstat(logit.mod5,resp=hsb$program)
```

#### For ROC w/ multivariate logitistic regression
```{r roc5}
roc.data5 <-roc(logit.mod5$y , logit.mod5$fitted.values,ci=T)
names(roc.data5)
plot((1-roc.data5$specificities),roc.data5$sensitivities,
      type='l',
	  lwd=2,
	  col="blue",
     main="ROC for Models 5, 1 and chance",
	 xlab="1-sensitivity (false positives)",
	 ylab="Specificity (true positives)"
	 )
roc.data1 <- roc(logit.mod1$y, logit.mod1$fitted.values,ci=T)
lines((1-roc.data1$specificities), roc.data1$sensitivities, col="pink", lwd=2)
text(.1,.99,"model 5: c=.79",col="blue")
text(.1,.94,"model 1: c=.75",col="pink")
text(.1,.89,"chance:  c=.50",col="black")
lines(c(0,1),c(0,1), lty=2)
legend("bottomright",
        legend=c("math + ses + schyp + ses*schtyp", "math","chance"),
		col=c("blue","pink","black"), lty=c(1,1,2), lwd=2)
```

# More predictors and problem of multicolinearity                  

```{r model_6}
summary(logit.mod6  <- glm(program ~ math + ses + sctyp + ses*sctyp + ses*math + math*sctyp + math*ses*sctyp, data=hsb,family = binomial)) 
```

## center all variables
```{r cfat_model}
hsb$cmath <- scale(hsb$math,center=TRUE,scale=FALSE)
hsb$cses <- scale(hsb$ses,center=TRUE,scale=FALSE)
hsb$sctyp <- as.numeric(hsb$sctyp)
hsb$csctyp <- scale(hsb$sctyp,center=TRUE,scale=FALSE)

summary(logit.mod7 <- glm(program ~ cmath + cses + csctyp + cses*csctyp + cses*cmath + cmath*csctyp +cmath*cses*csctyp, data=hsb,family = binomial)) 

summary(logit.mod8 <-  glm(program ~ math + ses + sctyp + ses*sctyp + rdg + sci + wrtg + civ, data=hsb,family = binomial)) 
```

## Examine correlations among achievement tests
```{r corrs}
var.names <- c("math","rdg","sci","wrtg","civ")
var.data <- hsb[var.names]
cor.mat <- cor(var.data,method="pearson")
round(cor.mat,2)

# Note eigen-vectors can be reflected (i.e., multiply by -1)
eigen(cor.mat)

hoslem.test(logit.mod5$y, fitted(logit.mod5), g=10)
```
## qq plot residuals

Using Pearson (adjusted standardized residuals
```{r}
resid <- rstandard(logit.mod5, pearson=TRUE)
qqnorm(resid,
     main="QQ plot of  Pearson  Residuals",
	   col="blue")
lines(resid,resid)
```

## ROC and concordance#

Concordance indices
```{r}
Cstat(logit.mod1,resp=hsb$program)
Cstat(logit.mod2,resp=hsb$program)
Cstat(logit.mod3,resp=hsb$program)
Cstat(logit.mod4,resp=hsb$program)
Cstat(logit.mod5,resp=hsb$program)
```

## For ROC w/ multivariate logitistic regression
```{r roc5b}
roc.data5 <-roc(logit.mod5$y , logit.mod5$fitted.values,ci=T)
names(roc.data5)
plot((1-roc.data5$specificities),roc.data5$sensitivities,
      type='l',
	  lwd=2,
	  col="blue",
     main="ROC for Models 5, 1 and chance",
	 xlab="1-sensitivity (false positives)",
	 ylab="Specificity (true positives)"
	 )
roc.data1 <- roc(logit.mod1$y, logit.mod1$fitted.values,ci=T)
lines((1-roc.data1$specificities), roc.data1$sensitivities, col="pink", lwd=2)
text(.1,.99,"model 5: c=.79",col="blue")
text(.1,.94,"model 1: c=.75",col="pink")
text(.1,.89,"chance:  c=.50",col="black")
lines(c(0,1),c(0,1), lty=2)
legend("bottomright",
        legend=c("math + ses + schyp + ses*schtyp", "math","chance"),
		col=c("blue","pink","black"), lty=c(1,1,2), lwd=2)
```



# Even more More predictors and problem of multicolinearity                    #

```{r}
summary(logit.mod6  <- glm(program ~ math + ses + sctyp + ses*sctyp + ses*math + math*sctyp + math*ses*sctyp, data=hsb,family = binomial)) 

# center all variables
hsb$cmath <- scale(hsb$math,center=TRUE,scale=FALSE)
hsb$cses <- scale(hsb$ses,center=TRUE,scale=FALSE)
hsb$sctyp <- as.numeric(hsb$sctyp)
hsb$csctyp <- scale(hsb$sctyp,center=TRUE,scale=FALSE)

summary(logit.mod7 <- glm(program ~ cmath + cses + csctyp + cses*csctyp + cses*cmath + cmath*csctyp +cmath*cses*csctyp, data=hsb,family = binomial)) 

hsb$sctyp <- as.factor(hsb$sctyp)
summary(logit.mod8 <-  glm(program ~ math + ses + sctyp + ses*sctyp + rdg + sci + wrtg + civ, data=hsb,family = binomial)) 

screenreg(list(logit.mod1,logit.mod5, logit.mod8),
          custom.model.names=c("Model 1", "Model 5", "Model 8"))

# Examine correlations among achievement tests
var.names <- c("math","rdg","sci","wrtg","civ")
var.data <- hsb[var.names]
cor.mat <- cor(var.data,method="pearson")
round(cor.mat,2)

# Note eigen-vectors can be reflected (i.e., multiply by -1)
eigen(cor.mat)
```
