R/ad.R

In andrewthomasjones/adtest:

#d^2 ij
d2ij<-function(y,s2){
  ((y-mean(y))^2)/s2
}

wstarij <- function(y,s2,nu){
  n=length(y)
  return((nu*n*d2ij(y,s2)/((nu+1)*(n-1)-n*d2ij(y,s2))))
}


#function to calculate aj
aij<-function(y,s2,nu){
  return(pf(wstarij(y,s2,nu),1,nu,lower.tail=F))
}


AD_main <- function(y) {

  ni<-dim(y)[1]
  g<-dim(y)[2]

  a<-array(0,dim(y))
  d<-array(0,dim(y))
  w<-array(0,dim(y))
  a_o<-array(0,dim(y))

  AD<-array(0,g)
  AD_1<-array(0,g)
  AD_2<-array(0,g)


  n<-ni*g #assume equal size groups

  nu=n-g-1

  #pooled variance
  s2=sum((ni-1)*apply(y,2,var)/(n-g))

  for(i in 1:g){
    d[,i]<-d2ij(y[,i],s2)
    w[,i]<-wstarij(y[,i],s2,nu)
    a[,i]<-aij(y[,i],s2,nu)
    a_o[,i]<-a[order(a[,i]),i]

    AD[i]<-(-ni)

    for(j in 1:ni){
      AD_1[i]<- AD_1[i] + -(3/ni)^(1/2)*(2*a_o[j,i]-1)
      AD_2[i]<- AD_2[i] + -(5/ni)^(1/2)*(1/2)*(3*(2*a_o[j,i]-1)^2-1)
      AD[i]<- AD[i] -(1/ni)*(2*j-1)*(log(a_o[j,i])+log(1-a_o[(ni-j+1),i]))
    }

  }

  return(list(AD=AD,AD_1=AD_1, AD_2=AD_2))
}



ad.test<-function(y, pooled=F){
  y<-as.matrix(y)
  g<-dim(y)[2]
  out<-NA

  #For each value of g, it outputs . A, A_1, and A_2 (first two components of A) ,
  #and the value of A and A-1,A-2 for the areas combined if g>1

  if((g>1 & pooled==F)| g==1){

    AD_sep<-array(0,g)
    pooled=F

    for(i in 1:g){
      AD_sep[i]<-AD_main(as.matrix(y[,i], ncol=1))[1]
    }

    out<-matrix(AD_sep, ncol=1)

    if(!is.null(colnames(y))){
      row.names(out)<-paste0("Group_", (colnames(y)))
    }else{
      row.names(out)<-paste0("Group_", (1:g))
    }

    colnames(out)<-c("A")

  }


  if(g>1 & pooled==T){

    AD_areas = AD_main(y)
    AD_combined=AD_main(matrix(y,ncol=1))
    out_list<-c(AD = AD_areas ,AD_combined=AD_combined)
    out<-as.data.frame(AD_areas)
    out<-rbind(out,AD_combined)
    colnames(out)<-c("A", "A_1", "A_2")

    if(!is.null(colnames(y))){
      row.names(out)<-c(paste0("Group_", (colnames(y))), "Combined")
    }else{
      row.names(out)<-c(paste0("Group_", (1:g)), "Combined")
    }

  }

  return(out)
}