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))