---
title: "Bully Example from Lecture"
author: "Carolyn J. Anderson"
date: '2022-10-05'
output: 
   html_document:
     toc: yes
---

This should reproduce the bully example results that are given in lecture.  The results in the lecture notes were created using SAS.  Some of the results below (i.e., marginal distribution graphs) seem to be a bit off.  

```{r setup}
library(MASS)  # for negative binomial
library(pscl)  # for zero inflated

setwd("D:/Dropbox/edps 589/4 glm")

bully <- read.table("bully_data.txt",header=TRUE)
```
# Exploratory

### Some basic descriptive statistics
```{r}
(n <- length(bully$nom))
(ybar <- mean(bully$nom))
(s <- sd(bully$nom))
(s2 <- var(bully$nom))
summary(bully$nom)
```

### Histogram of bully nominations

```{r}
hist(bully$nom,
     main="Distribution of Number of Bully Nominations",
     xlab="Number of Nominations that a student is Bully",
     breaks=50,
     col="cyan")
text(41.2,200,"N = 291")	 
text(40,190,"Mean = 2.49")	 
text(41,180,"sd = 6.42")
text(41.1,170,"s2 = 41.26")
text(39.1,160,"min = 0")
text(39.1,150,"max = 50")
```

### Histogram of bully nominations with null poisson overlay
```{r}
y <- seq(1:51)-1
d.poisson <- dpois(y, ybar)	 
	 
hist(bully$nom, prob=TRUE, col="cyan", breaks=50, 
     main="Bully Nominations w/ Poisson Overlayed",
     xlab="Number of Bully Nominations")     
#prob=TRUE for probabilities not counts
lines(y,d.poisson, col="darkgreen", lty=2,  lwd=1.5)   # added a density estimate 
```

# Poisson loglinear regression 
Although we don't expect this to fit...

```{r}
mod1 <- glm(nom ~ bsc, data=bully, family=poisson)
summary(mod1)
```

Let's see how well this is doing
```{r}
bully$pred1 <- fitted.values(mod1)
pred <- predict(mod1, newdata = NULL,
            type = c("link", "response", "terms"),
            se.fit = TRUE
        )
bully$lower1 <- exp(pred$fit  - 1.96*pred$se.fit)
bully$upper1 <- exp(pred$fit  + 1.96*pred$se.fit)
bully$pred1  <- exp(pred$fit)

bully.sort <- bully[order(bully$bsc),]
plot(bully.sort$bsc,bully.sort$nom,type='p',pch=20, col="blue",
     main='Data, Fitted Line, and 95% Confidence', 
     xlab='Bully Scale Score', ylab='Fitted Count')
lines(bully.sort$bsc,bully.sort$pred1,col="darkgreen",lty=1,lwd=2)
lines(bully.sort$bsc,bully.sort$lower1,col="darkgreen",lty=2)
lines(bully.sort$bsc,bully.sort$upper1,col="darkgreen",lty=2)
```


# Other models   -- digression  
### identity link in poisson
```{r}
summary(mod.Pidentity <- glm(nom ~ bsc, data=bully, family=poisson("identity")))
mod.Pidentity$aic
```

### identity link and negative binomial
```{r}
summary(nb.indentity <- glm.nb(nom ~ bsc, data=bully, link=identity))
nb.indentity$aic
```
But the line is for fitted mean..so another way to look at the goodness of fit is to group data based on bully scale score and re-fit the model, etc.

This is ONLY for assessing fit (fit to individual level is what should be report and interpreted

Group data and get a look-see at how good it may be fitting the data

```{r}
bully$grp <- cut(bully$bsc, quantile(bully$bsc, c(0, 1/7, 2/7, 3/7, 4/7, 5/7, 6/7, 1))) 	
grp.bsc <- aggregate(bully$bsc, by= list(bgrp=bully$grp), FUN=mean)

grp.nom <- aggregate(bully$nom, by= list(bgrp=bully$grp), FUN=mean)

grpdata <- cbind(grp.bsc, grp.nom)
grpdata <- as.data.frame(grpdata)
names(grpdata) <- c("grp","mbsc","grp2","mnom")
plot(grpdata$mbsc,grpdata$mnom)
```

```{r}
mod1b <- glm(mnom ~ mbsc, data=grpdata, family=poisson)
summary(mod1b)
#-- test of fit
1 - pchisq(mod1b$deviance,mod1b$df.residual)   # look good

#  S3 method for class 'glm'
pred <- predict(mod1b, newdata = NULL,
            type = c("link", "response", "terms"),
            se.fit = TRUE
        )
```
This looks difference from SAS but numerical results are spot on (i.e. equal for all quanties so ...

```{r}
lower <- exp(pred$fit  - 1.96*pred$se.fit)
upper <- exp(pred$fit  + 1.96*pred$se.fit)
fit   <- exp(pred$fit)
```

```{r}
plot(grpdata$mbsc,grpdata$mnom, 
     pch=20, col="blue", 
     main='Grouped: Means bully scale by Nominations (95% Confidence)',
     xlab='mean bully scale score',
     ylab='observed and fitted number nominations',
     ylim=c(0,35),
     xlim=c(1,4) )
lines(grpdata$mbsc,fit,col='darkgreen',lty=1,lwd=2)
lines(grpdata$mbsc,lower,col='darkgreen',lty=2)
lines(grpdata$mbsc,upper,col='darkgreen',lty=2)
```



Fitted Marginal distribution (back to individual level)
```{r}
ss <- seq(1:length(bully$bsc))
yvalues <- seq(1:51)
fitprob <- matrix(99,nrow=291,ncol=51)
for (i in ss) {
  for (yindex in yvalues) { 
    yobs <- yindex-1
    mu <- exp(-0.66308 + 0.81434*bully$bsc[i] )
    fitprob[i,yindex] = ( exp(-mu) * mu**yobs)/(factorial(yobs))
  }
}
# summing over y values (colums) should all equal 1
rowSums(fitprob)

# mean over kids (rows) is what I want to plot
marginals <- colSums(fitprob)/length(bully$bsc)

# Histogram of bully nominations with null poisson overlay
y <- seq(1:51)-1
d.poisson <- dpois(y, ybar)	 
	 
hist(bully$nom, prob=TRUE, col="cyan", breaks=50, 
     main="Bully Nominations with Poisson Overlayed",
     xlab="Number of Bully Nominations")     # prob=TRUE for probabilities not counts
lines(y,d.poisson, col="darkgreen", lty=2,  lwd=1)     # added a density estimate 
lines(y,marginals, col="black", lty=1, lwd=2)          # mod 1, poisson with predictor
text(30,0.68,"Data",col="cyan")
text(36,0.62,"Poisson no predictors",col="darkgreen")
text(39.4,0.66,"Poisson nom ~ bullyscale score",col="black")
```

# Negative Binomial

```{r}
library(MASS)
```

## Fit model to individual level data
```{r}
mod2 <- glm.nb(nom ~ bsc, data=bully)  # dispersion phi = 1/"disperson from R"
summary(mod2)
```

Using my grouped data
```{r}
mod2b <- glm.nb(mnom ~ mbsc, data=grpdata)
summary(mod2b)
#-- test of fit
1 - pchisq(mod2b$deviance,mod2b$df.residual)   
# looks good
```

##  S3 method for class 'glm'
```{r}
pred <- predict(mod2b, newdata = NULL,
            type = c("link", "response", "terms"),
            se.fit = TRUE         )

lower <- exp(pred$fit  - 1.96*pred$se.fit)
upper <- exp(pred$fit  + 1.96*pred$se.fit)
fit   <- exp(pred$fit)
```

And a graph
```{r}
plot(grpdata$mbsc,grpdata$mnom, 
     pch=20, col="blue", 
     main='NegBin to Grouped: Means bully scale by Nominations (95%)',
     xlab='mean bully scale score',
     ylab='observed and fitted number nominations',
     ylim=c(0,35),
     xlim=c(1,4) )
lines(grpdata$mbsc,fit,col='darkgreen',lty=1,lwd=2)
lines(grpdata$mbsc,lower,col='darkgreen',lty=2)
lines(grpdata$mbsc,upper,col='darkgreen',lty=2)                
```
Doesn't look same as SAS, but not sure why.


## Fitted Marginal distribution 

Using back to individual level 
 
This seems to be off a bit off; not sure why.

```{r}
ss <- seq(1:length(bully$bsc))
yvalues <- seq(1:51)
fitprob2 <- matrix(99,nrow=291,ncol=51)               # For negative binomial
phi = 0.2858
for (i in ss) {
  for (yindex in yvalues) { 
    yobs <- yindex-1
    mu <- exp(-1.1940 + 1.0921*bully$bsc[i] )
    term1 <- lgamma(mu+phi)/(factorial(yobs)*lgamma(phi))
    term2 <- ( phi/(mu+phi))**phi 
    term3 <- ( mu/(mu+phi))**yobs  
    fitprob2[i,yindex] = term1*term2*term3
  }
}
```
 
Mean over kids (rows) is what I want to plot
```{r}
marginals2 <- colSums(fitprob2)/length(bully$bsc)

# Histogram of bully nominations with null poisson overlay
y <- seq(1:51)-1
d.poisson <- dpois(y, ybar)	 
	 
hist(bully$nom, prob=TRUE, col="cyan", breaks=30, 
     main="Bully Nominations with Fitted Overlayed",
     xlab="Number of Bully Nominations",
     ylim=c(0,1)
     )     # prob=TRUE for probabilities not counts
lines(y,d.poisson, col="darkgreen", lty=2,  lwd=1)         # add a density estimate 
lines(y,marginals, col="black", lty=1, lwd=1)              # mod 1, poisson with predictor
lines(y,marginals2,col="red",lty=1, lwd=2)                 # negative binomial
text(30,0.30,"Data",col="cyan")                            # zip
text(31,0.25,"Poisson no predcitors",col="darkgreen")
text(32.4,0.20,"Poisson nom ~ bullyscale score",col="black")
text(32.4,0.15,"NegBin nom ~ bullyscale score",col="red")
```

These look a bit off.

# ZIP #######################

```{r}
library(pscl)
```

```{r}
summary(mod3 <- zeroinfl(nom ~ bsc | 1, data = bully))
```

```{r}
summary(mod4 <- zeroinfl(nom ~ bsc | bsc, data = bully))
```

```{r}
ss <- seq(1:length(bully$bsc))
yvalues <- seq(1:51)
fitprob4 <- matrix(99,nrow=291,ncol=51)               # 
# For zip (mod4)

for (i in ss) {
  for (yindex in yvalues) { 
    yobs <- yindex-1
    mu.zip <- exp(0.4967 + 0.59609*bully$bsc[i] )
    pie <- exp(1.4967 - 0.7716*bully$bsc[i])/( 1 + exp(1.4967 - 0.7716*bully$bsc[i]))
    if (yobs==0) {
       fitprob4[i,yindex] <- pie + (1-pie)*exp(-mu.zip)
    } else {
       fitprob4[i,yindex] <- (1-pie)*exp(-mu.zip)*(mu.zip**yobs)/factorial(yobs)
    }
  }
}
```

Mean over kids (rows) is what I want to plot

```{r}
marginals4 <- colSums(fitprob4)/length(bully$bsc)
y <- seq(1:51)-1
d.poisson <- dpois(y, ybar)	 
	 
hist(bully$nom, prob=TRUE, col="cyan", breaks=30, 
     main="Bully Nominations with Fitted Overlayed",
     xlab="Number of Bully Nominations",
     ylim=c(0,1)
     )     # prob=TRUE for probabilities not counts
lines(y,d.poisson, col="darkgreen", lty=2,  lwd=1)         # add a density estimate 
lines(y,marginals, col="black", lty=1, lwd=1)              # mod 1, poisson with predictor
lines(y,marginals2,col="red",lty=1, lwd=2)                 # negative binomial
lines(y,marginals4,col="violet", lty=4, lwd=2)
    # zip
text(30,0.35,"Data",col="cyan")  
text(32.4,0.30,"Poisson nom ~ bullyscale score",col="black")
text(32.4,0.25,"NegBin nom ~ bullyscale score",col="red")
text(32.4,0.20,"ZIP nom ~ bullyscale score",col="violet")
```

The marginal predictions for zip and negbin seem to be a bit off.  