### Week 8: Count Dependent Variables
### Code from class

require(foreign)

## classical OLS versus MLE linear-normal
dail <- read.dta("dailcorrected.dta")
summary(m1.ols <- lm(votes1st ~ spend_total*incumb, data=dail))
summary(m1.mle <- glm(votes1st ~ spend_total*incumb, data=dail))
summary(m1.ols)$sigma^2
require(lmtest)
lrtest(m1.mle)

# trying to compute AIC
m1.mle$aic
(n2ll <- -2*logLik(m1.mle))
2*m1.mle$rank + n2ll

# likelihood ratio test for bigger model (and compare to F-test)
summary(m2.mle <- glm(votes1st ~ spend_total*incumb + electorate + senator,
                      data=dail))
require(lmtest)
lrtest(m1.mle, m2.mle)
summary(m2.ols <- lm(votes1st ~ spend_total*incumb + electorate + senator,
                     data=dail))
anova(m1.ols, m2.ols)


## illustrate shape of exponential curves
par(mar=c(4,4,1,1))
x <- matrix(cbind(1,seq(-1,5,.1)), ncol=2)
plot(x[,2], exp(x%*%c(2,-.5)), type="l", xlab="X", ylab="exp(XB)")
text(.3,9,"b1=2, b2=-.5")
lines(x[,2], exp(x%*%c(0,.7)), type="l")
text(2.5,9.3,"b1=0, b2=.7")
lines(x[,2], exp(x%*%c(-2,1)), type="l")
text(4.3,5,"b1=-1, b2=.1")


## illustrate differences of Poisson and negbin using Benoit (1996)
## From Table 2 weede <- read.dta("weede.dta") weede <-
weede <- read.dta("weede.dta")
weede <- weede[complete.cases(weede$ssal6080, weede$poldem65, weede$lpopln70,
                              weede$lmilwp70, weede$ecintdep, weede$lecnpc70),]
summary(rdwi.ols <- lm(ssal6080 ~ poldem65, data=weede))
summary(rdwi.pois <- glm(ssal6080 ~ poldem65, family=poisson,
                         data=weede))
2*logLik(rdwi.pois)
require(MASS)
summary(rdwi.nb <- glm.nb(ssal6080 ~ poldem65, data=weede))

rdwi.pois2 <- glm(ssal6080 ~ poldem65+lecnpc70+ecintdep,
                  family=poisson, data=weede)
lrtest(rdwi.pois, rdwi.pois2)

# model 1 from table 4
summary(rdwi.pois.m1tab4 <-
        glm(ssal6080 ~ poldem65+lpopln70+lmilwp70+ecintdep+lecnpc70,
            family=poisson, data=weede))
# using negative binomial model
summary(rdwi.pois.m1tab4 <- glm.nb(ssal6080 ~
        poldem65+lpopln70+lmilwp70+ecintdep+lecnpc70, data=weede))


## Example of logit using odds-ratios as coefficients
require(foreign)
dail <- read.dta("dailcorrected.dta")
won.logit <- glm(wonseat ~ incumb*spend_total + electorate,
                 family=binomial(link="logit"), data=dail)
summary(won.logit)$coefficients
exp(summary(won.logit)$coefficients)
exp(confint(won.logit))

## in Stata:
# use dailcorrected, clear
# logit wonseat incumb spend_total spendXinc electorate, or nolog


## Fitted values
summary(p1 <- glm(ssal6080 ~
                  fh73+lpopln70+lmilwp70, family=poisson, data=weede))
x <- seq(2,14,1)
bhat <- p1$coefficients
xfitted <- cbind(rep(1,length(x)), x,
                 mean(weede$lpopln70), mean(weede$lmilwp70))
xbfitted <- xfitted %*% matrix(bhat,nrow=length(bhat))
yfitted <- exp(xbfitted)
plot(x, yfitted, type="b")
cbind(x,yfitted)