## AIC & BIC(for monotone missing data & variable order p)

PL<-function(data, covariate, BIC){
  data <- as.matrix(data)
  Z <- as.matrix(covariate)
  n <- ncol(data)
  p <- rep(0, n)
  
  for(i in 2:n){

    Z <- as.matrix(Z[!is.na(data[,i]),])    
    data <- data[!is.na(data[,i]),]      #delete missing data
    N <- nrow(data)  #Ni

    if(BIC==TRUE)  c.N <- log(N) 
    else  c.N <-2

    p1 <- 0
    for(j in 1:(i-1)){
      A <- data[, i]

      if(p1==0)    
        B <- Z
      else
        B <- cbind(Z, as.matrix(data[, (i-min(p1, i-1)):(i-1)]))     

      C <- cbind(Z, as.matrix(data[, (i-min(j, i-1)):(i-1)]))     
 
      RSS1 <- log(var(residuals(lm(A~B)))*(N-1))+p1*c.N/N
      RSS2 <- log(var(residuals(lm(A~C)))*(N-1))+j*c.N/N

      if(RSS1 > RSS2)  p1 <- j
    }
    p[i] <- p1
  }
  p
}

## Speech recogintion data
speech<-read.table("cochlear09.dt")[,2:5]
covariate<-rbind(t(matrix(rep(c(1,0),20), nrow=2)), t( matrix(rep(c(0,1),21), nrow=2)))
> PL(speech, covariate, BIC=TRUE)
[1] 0 1 1 1
> PL(speech, covariate, BIC=FALSE) ##AIC
[1] 0 1 1 1



## Treatment A cattle growth data
cattle<-read.table("kenwarda.dt")
covariate<-matrix(rep(1, 30), ncol=1)
> PL(cattle, covariate, BIC=TRUE)
 [1] 0 1 1 1 1 1 1 2 2 1 1
> PL(cattle, covariate, BIC=FALSE) ##AIC
 [1] 0 1 1 1 1 1 1 2 3 1 1



## 100-km race data
race100k<-read.table("race100k.dt")  ##complete set of 80 competitors 
race<-race100k[, 2:11]
covariate<-matrix(rep(1, 80), ncol=1)
> PL(race, covariate, BIC=TRUE)
 [1] 0 1 1 1 2 1 1 2 3 5
> PL(race, covariate, BIC=FALSE)
 [1] 0 1 1 1 3 1 2 2 3 5