#set of functions to accompany paper by Sin-Ho Jung1 and Chul W. Ahn
#Sample size for a two-group comparison of repeated binary
#measurements using GEE
#STATISTICS IN MEDICINE
#Statist. Med. 2005; 24:2583–2596

#logit link
logit<-function(p){
   if (abs(p-.5) >= .5) {
      stop("probability must be between 0 and 1")
   }
   log(p)-log(1-p)
}

#inverse logit link
logiti<-function(a, b=0, t=0) exp(a + b*t)/(1 + exp(a + b*t))

#absulute difference between 
#i'th element of the list and 
#every other element in the list
absdiff<-function(i,lst=c(0)){
   len<-length(lst)
   if (i < 1) {
      stop("index must be above 0")
   }
   if (i > len) {
      stop("index must be not greater than length of the list")
   }
   abs(lst-lst[i])
}

#create AR(1) structured correlation matrix
ar1corr<-function(rho,dim){
   if (abs(rho) >= 1) {
      stop("Parameter in AR(1) structure must be between -1 and 1")
   }
   if (dim < 1) {
      stop("dimension of the matrix must be positive integer")
   }
   matr<-rho^sapply(1:dim,absdiff,lst=1:dim)
   matr
}
#create compound symmetry correlation matrix
cscorr<-function(rho,dim){
   if (abs(rho) >= 1) {
      stop("Parameter in AR(1) structure must be between -1 and 1")
   }
   if (dim < 1) {
      stop("dimension of the matrix must be positive integer")
   }
   matr<-matrix(rep(rho,dim*dim),nrow=dim)+as.matrix(diag(1-rho,dim))
   matr
}

#independent missing structure
# delta_ij=delta_i*delta_j
indMCAR<-function(obs)
   obs%*%t(obs)

#pairwise minimum value between
#i'th element of the list and 
#every other element in the list
minim<-function(i,lst=c(0)){
   len<-length(lst)
   if (i < 1) {
      stop("index must be above 0")
   }
   if (i > len) {
      stop("index must be not greater than length of the list")
   }

   pmin(lst,rep(lst[i],6))
}

#monotone missing
# delta_ij=delta_min(i,j)
monMCAR<-function(obs){
   len<-length(obs)
   sapply(1:len,minim,lst=obs)
}

#at this moment corStruct="cs" or "ar1" 
#or "ud" - user defined, must provide correlation matrix
#at this moment miss="ind" or "mon" - independent or monotone
sampleSize<-function(ac,bc,bt,t,obs,rho,
          corStruct="ar1",miss="ind",at=ac,alpha=.05,beta=.2,r1=.5){
   if(abs(r1-.5)>.5){
      stop("Proportion for the group 1 should be strictly between 0 and 1.")
   }else{
      r2<-1-r1
   }
   if(length(t)==1)t<-0:t
   if(length(obs)==1)obs<-obs^t
   if(length(t)!=length(obs))
      stop("Percent of observed values must be assigned for each time point.")
   d<-abs(bc-bt)
   tr <- logiti(at, bt, t)
   co <- logiti(ac, bc, t)
   if(corStruct=="cs"){
      rho<-cscorr(rho,length(t))
   }else if(corStruct=="ar1"){
      rho<-ar1corr(rho,length(t))
   }else if(corStruct=="ud"){
      if((nrow(rho)!=length(obs))||(ncol(rho)!=length(obs)))
         stop("Need a square matrix with dimensions equal to number of time points")
   }else{
      stop("Such correlation structure is not implemented")
   }
   if(miss=="ind"){
      obsm<-indMCAR(obs)
   }else if(corStruct=="ar1"){
      obsm<-monMCAR(obs)
   }else{
      stop("Such missing scheme is not implemented")
   }
   tau1<-sum(obs*co*(1-co)*t)/sum(obs*co*(1-co))
   s12<-sum(obs*co*(1-co)*(t-tau1)^2)
   c12prep<- obsm*rho*sqrt((co*(1-co))%*%t(co*(1-co)))*((t-tau1)%*%t(t-tau1))
   c12<-sum(c12prep-diag(diag(c12prep)))
   v1<-(c12+s12)/s12^2
   tau2<-sum(obs*tr*(1-tr)*t)/sum(obs*tr*(1-tr))
   s22<-sum(obs*tr*(1-tr)*(t-tau2)^2)
   c22prep<- obsm*rho*sqrt((tr*(1-tr))%*%t(tr*(1-tr)))*((t-tau2)%*%t(t-tau2))   
   c22<-sum(c22prep-diag(diag(c22prep)))
   v2<-(c22+s22)/s22^2
   n<-(qnorm(alpha/2)+qnorm(beta))^2*(v1/r1+v2/r2)/d^2
   frac1<-1/(1+sqrt(v2/v1))
   list(n=ceiling(n),frac1=frac1)
}

#---test sample size calc
ac<-logit(.75)
bc<- (logit(.5)-logit(.75))/5
bt<-0
t<-0:5
obs<- c(1, 0.95, 0.9, 0.85, 0.8, 0.75)
rho<-.8
res1ind<-sampleSize(ac,bc,bt,t,obs,rho)
res1ind
res2ind<-sampleSize(ac,bc,bt,t,obs,rho,r1=res1ind$frac1)
res2ind
res3mon<-sampleSize(ac,bc,bt,t,obs,rho,miss="mon")
res3mon
res4mon<-sampleSize(ac,bc,bt,t,obs,rho,miss="mon",r1=res3mon$frac1)
res4mon




#tests
logit(.75)
logiti(0)
logiti(1)
logiti(1,2,1)
t<-0:5
absdiff(3,t)
ar1corr(.8,6)
ar1corr(-.8,6)
cscorr(.8,6)
minim(3,t)
obs<- c(1, 0.95, 0.9, 0.85, 0.8, 0.75)
indMCAR(obs)
monMCAR(obs)