#open R and installing the "geepack" package
install.packages("geepack")
library(geepack)
#import data from text file
data<-read.table(file.choose(),header=T,fill=T)
attach(data)
#explore the variable name
names(data)
#summary stat
summary(data)
#finding the variation of variable(s) > aaa= variable name of each (continous data)
sd(aaa)
#calculating multivariable analysis using generalised estimating eqautions (GEE) model for binary longitudinal data
#univariate >> run model using independent variable for each then
data.gee<-geeglm(dis~f1,data=data,id=pondcode,family=binomial("logit"),corstr="exchangeable")
#multivariate >> run model until all remaining factors were significant
#easily remove non-selected varibale by untype the variable(s) name in the GEE model. The removal method should remove most non-significant varaible firstly 
data.gee<-geeglm(dis~f5+f7+f10+f16+f19+f53+f75+f77+f87+f89+f93,data=data,id=pondcode,family=binomial("logit"),corstr="exchangeable")
summary(data.gee)

summary(data.gee)
##finding odd ratio from coefficient and its 95%CI
exp(estimatevalue)
exp((estimatevalue)+(1.96*robustse*estimatevalue));exp((estimatevalue)-(1.96*robustse*estimatevalue))
#compare with alpha 0.05 using critical z (1.96) >> robust z > 1.96 or < -1.96 is significant
#calculating 95%CI of estimate
se<-summary(data.gee)$coefficients["variable name","Robust S.E."]
coef(data.gee)["vairbale name"]+c(-1,1)*se*qnorm(0.975)
#model comparison > comparing the final model with naive model
data.gee2<-geeglm(disstatus~xxx1+xxx2+...+xxxn,data=data,id=pondcode,family=binomial("logit"),corstr="")
##creating the graph using "ggpubr" package > selective variable(s)
install.packages("")
library()

#model vlidation using "geepack" packages >> QIC >> QIC value of exchangeable should be lower than another one
data.gee.qic<-geeglm(disstatus~xxx1+xxx2+...+xxxn,data=data,id=pondcode,family=binomial("logit"),corstr="exchangeable")
data.gee.qic2<-geeglm(disstatus~xxx1+xxx2+...+xxxn,data=data,id=pondcode,family=binomial("logit"),corstr="independence")
QIC(data.gee.qic);QIC(data.gee.qic2)