# homework 4
3 Edps 589
# Fall 2018
# C.J.Anderson


setwd("C:\\Users\\cja\\Dropbox\\edps 589\\homework")

# problem 1:  alchol consumption and malformation

a <- read.table("hm4_problem1_data.txt",header=TRUE)
a$total <- a$absent + a$present;

# (a) Linear probability model
summary(lin.mod <- glm(present/total ~ score, family=binomial("identity"),data=a) )


# (b) probit model
summary(probit.mod <- glm(present/total ~ score, family=binomial("probit"),data=a) )

# (c) logit model
summary(logit.mod <- glm(present/total ~ score, family=binomial("logit"),data=a) )

a$p <- a$present/a$total

plot(a$score,a$p, type="p",col="black", pch=20,
     main="Fitted Value: Infant Malformation and Alcohol Consumption",
	 xlab="Alcohol Consumption Score",
	 ylab="Malformation Present")
lines(a$score,logit.mod$fitted,col="red")
lines(a$score,probit.mod$fitted,col="blue")
lines(a$score,lin.mod$fitted,col="green")
legend(0.0,0.026, legend=c("Observed Proporitions", "Logit", "Probit","Linear probability"),
       col=c("black","red", "blue","green"), lty=1, cex=0.8)
	   
###################################################
# Problem 2                                       #
# From Agresti, 3.11                              #
###################################################
# Parts (a) and (b)

# Read in data
(wafer <- read.table("a3_11_12_data.txt",header=TRUE))

# Fit model using treatment as predictor
summary( poi.mod1 <- glm(imperf ~ treat, data=wafer, family=poisson ) )

# Getting linear predictors
(log.muB <- poi.mod1$coefficients[1] + poi.mod1$coefficients[2])
(log.muA <- poi.mod1$coefficients[1] )

# Estimate of beta hat
(beta.hat <- log.muB - log.muA)

# For interpretation
(odds.ratio <- exp(beta.hat))


# Part (c)
# Wald statistic -- square z
(wald <- 3.332**2)
1-pchisq(wald,1)


# Likelihood ratio test
(lr <- poi.mod1$null.deviance-poi.mod1$deviance)
1-pchisq(lr,1)

# Part (d)
# -- 99% for beta
(z99 <- qnorm(.995))
(lower <- poi.mod1$coefficients[2] - z99*.1764)
(upper <- poi.mod1$coefficients[2] + z99*.1764)

# -- 99% for muB/muA
(mu.lo <- exp(lower))
(mu.hi <- exp(upper))

##### Alternate answer ways to answer problem 2 ####
(table(wafer$treat))

(p <- 50/140)

# large sample z
(z <- (p - .5)/sqrt(.5*.5/140) )

# Or exact
binom.test(50,140,p=.5,conf.level=.99)


###################################################
# Problem 3                                       #
# 3.12 Agresti                                    #
###################################################
summary(poi.mod2 <- glm(imperf ~ treat + thick, data=wafer, family=poisson ) )
anova(poi.mod2)

1-pchisq(16.268,18)