options(contrasts=c("contr.treatment", "contr.poly"))

ly.sat <- glm(y ~ Cell * Remis * Sex, data=lymphoma, family=poisson)
summary(ly.sat)


ly.glm1 <- update(ly.sat, .~. - Cell:Remis:Sex)
summary(ly.glm1)
anova(ly.glm1, test="Chisq")



ly.glm2 <- update(ly.glm1, .~. - Remis:Sex)

summary(ly.glm2)

?predict.glm

pcount <- predict(ly.glm2, type="response")
lymphoma$pcount <- pcount # See the lymphoma dataset
data.frame(observed=lymphoma$y, fitted=pcount)
fit.prob <- pcount/(sum(pcount)) 
lymphoma$fitprob <- fit.prob # See the lymphoma dataset
fit.prob



ly.pois <- glm(y ~ Cell+ Sex+ Remis+ Cell:Remis + Cell:Sex + Sex:Remis, data=lymphoma, family=poisson)
summary(ly.pois)


makepropdata  <- function() { 
     Cell <- c("nodular", "nodular", "diffuse", "diffuse")
     Sex <- c("male", "female", "male", "female")
     y <- c(4, 6, 1, 1)
     n <- c(5, 8, 13, 4)
     data.frame(Cell=Cell, Sex=Sex, y=y, n=n)
}
newlymphoma <- makepropdata()

ly.bino <- glm(y/n ~ Sex+Cell, data=newlymphoma, family=binomial, weights=n)

summary(ly.bino)


convbindata <- function()
{
	Cell <- rep(lymphoma[, 2], lymphoma[, 1])
	Sex <- rep(lymphoma[, 3], lymphoma[, 1])
	Remis <- rep(lymphoma[, 4], lymphoma[, 1])
	data.frame(Cell, Sex, Remis )
}
nlmph <- convbindata()
ly.bino2 <- glm(Remis~Sex+Cell, data=nlmph, family=binomial)
summary(ly.bino2)

coef(ly.pois)
coef(ly.bino)
coef(ly.bino2)


createdat <- function() { 
 y <- c(1, 2, 12, 3, 4, 6, 1, 1)
 Cell <- c("nodular", "nodular", "diffuse", "diffuse", "nodular", "nodular", "diffuse", "diffuse") 
 Sex <- c("male", "female", "male", "female", "male", "female", "male", "female")
 Remis <- c(rep("no", 4), rep("yes", 4))
data.frame(y, Cell, Sex, Remis)
}
lymphoma <- createdat()


# This is for the heartattack example 

convtotable  <- function() { 
     k <- length(heartattack$y)
     tmpy <- c(heartattack$y,  heartattack$n - heartattack$y)
     U <- c( rep("yes", k), rep("no", k) )
     V <- rbind.data.frame(heartattack, heartattack)
     V$U <- factor(U)
     V$Count <- tmpy 
     V
}

heart.pois <- glm(Count ~ Rtreat + Blocker + Site +   
              Rtreat:Blocker + Rtreat:Site + Blocker:Site + 
           U + U:Rtreat + U:Blocker + U:Site, data=newheart, family=poisson)
summary(heart.pois)

heart.bino <- glm(y/n ~ Rtreat+Blocker+Site,
     data=newheart, family=binomial, weights=n,subset=1:24)
summary(heart.bino)

# This was our hglm.final

coef(heart.pois)
coef(heart.bino)