rm(list=ls())

##########################
##### LAB. SESSION 2: ####
### FITTING LINEAR AND ###
### GENERALIZED LINEAR ###
######## MODELS ##########
##########################

# Slide 2

y~1 # Intercept only
y~x # Linear effect of x
y~x-1 # Removal of the intercept
y~x+I(x^2) # Quadratic effect of x

# Slide 3
state <- c(rep("Malawi",5),rep("Italy",2),rep("Madagascar",4))
state <- factor(state,levels=c("Italy", "Madagascar", "Malawi"))
contrasts(state)

# Slide 4

income <- c(rep(1,3),rep(2,5),rep(3,4))
income <- factor(income,levels=1:3)
contrasts(income)[lower.tri(contrasts(income))] <- 1
contrasts(income)

# Slide 5
beta0 <- -0.5
beta1 <- 1

sigma2 <- 0.5

n <- 100

set.seed(123)
x <- rnorm(n)
y <- beta0+beta1*x+rnorm(n,sd=sqrt(sigma2))
data.sim <- data.frame(y=y,x=x)

fit.lm <- lm(y~x,data=data.sim)

# Slide 6
plot(x,y)
abline(fit.lm)

z1 <- rnorm(n)
z2 <- rnorm(n)
data.sim2 <- data.frame(y=y,x=x,z1=z1,z2=z2)
fit.lm2 <- lm(y~x+z1+z2,data=data.sim2)

anova(fit.lm,fit.lm2)
F.test <- ((46.172-46.059)/2)/((46.059)/96)
F.test
p.value.F.test <- 1-pf(F.test,2,96)
p.value.F.test

# Slide 7
beta0 <- -0.5
beta1 <- 1

sigma2 <- 0.5

n <- 100

set.seed(123)
x <- rnorm(n)
eta <- beta0+beta1*x+rnorm(n,sd=sqrt(sigma2))
m <- rep(10,n)
y <- rbinom(n,size=m,prob=exp(eta)/(1+exp(eta)))
data.bin.sim <- data.frame(y=y,m=m,x=x)

fit.glm <- glm(cbind(y,m-y)~x,data=data.bin.sim,
           family=binomial(link="logit"))

# Slide 8
summary(fit.glm)

# Slide 9
z1 <- rnorm(n)
z2 <- rnorm(n)
data.bin.sim2 <- data.frame(y=y,m=m,x=x,z1=z1,z2=z2)
fit.glm2 <- glm(cbind(y,m-y)~x+z1+z2,
                data=data.bin.sim2,
                family=binomial(link="logit"))
logLik(fit.glm)
logLik(fit.glm2)
stat.test <- as.numeric(2*(logLik(fit.glm2)-logLik(fit.glm)))
p.value <- 1-pchisq(stat.test,2)
p.value

exp(coef(fit.glm2))

exp(confint(fit.glm2))