rm(list = ls(all = TRUE))
library(gee)
setwd("/home/jenny/practical_family/")   # change to your directory

# read plink ped file
ped  <-  read.table(file="rs2228145_plink.ped")
head(ped)
# add headers 
colnames(ped) <- c("FAMID", "PERSONID","FATH", "MOTH", "SEX","sIL6R","A1","A2")
head(ped)

# read plink covariate file
cov <- read.table(file="plink_covar.txt",header=T)
head(cov)
# add covariates to object ped
ped <- data.frame(ped,cov[,-c(1:2)])

# check the distribution
hist(ped$sIL6R,breaks=100)
# compute the mean and variance of sIL6R
mean(ped$sIL6R)
var(ped$sIL6R)


# to run the association test in gee, we first have to recode the genotypes to 0,1,2 (= the number of C alleles)
genonum <- rep(NA, nrow(ped))
genonum[which(ped$A1=="A" & ped$A2=="A") ] <- 0
genonum[which(ped$A1=="C" & ped$A2=="A") ] <- 1
genonum[which(ped$A1=="C" & ped$A2=="C") ] <- 2
# check the phenotype distribution for each genotype
boxplot(ped$sIL6R ~ genonum)

# run assocation test in gee
data  <- data.frame(ped,genonum)
# sort data by familynumber! (this is important)
data  <- data[order(data$FAMID),]
head(data)

# option independence

out <- gee(sIL6R~genonum + zage + PC1_NL + PC2_NL + PC3_NL + PC3_chip_effect + PC5_chip_effect + PC1_buccal,data=data, id=FAMID, family=gaussian, corstr="independence", maxiter=100, na.action=na.omit)
coeff <- summary(out)$coefficients
coeff

# option exchangeable

out <- gee(sIL6R~genonum + zage + PC1_NL + PC2_NL + PC3_NL + PC3_chip_effect + PC5_chip_effect + PC1_buccal,data=data, id=FAMID, family=gaussian, corstr="exchangeable", maxiter=100, na.action=na.omit)
coeff <- summary(out)$coefficients
coeff