########################################################################
#####                                                              #####
##### This R code runs the final model selected by QIC presented   #####
##### in Figures 3 and 4 and Table 3 of the paper                  #####
#####                                                              #####
##### ``Statistical Modeling of the Time Course of Tantrum Anger'' #####
##### by Peihua Qiu, Rong Yang, and Michael Potegal                #####
#####                                                              #####
##### Annals of Applied Statistics, accepted for publication in    #####
##### March 2009.                                                  #####
#####                                                              #####
########################################################################

####################################################################
outfile <- file("out.txt","w")   ### Write results to this file 

epsilon <- 10^-6                 ### A small number added to 
                                 ### diagonal elements of certain
                                 ### matrices to avoid singularity
####################################################################


####################################################################
dat <- matrix(scan("AOAS0801-028R1S4.dat"),ncol=10,byrow=T)
                                          ### read in the data

sub <- length(dat[,1])                    ### Number of observations
####################################################################


####################################################################
###                    Define certain functions
####################################################################
 
### Probability pi of a behavior in terms of all parameters

gabcd <- function(a0,a1,a2,b0,b1,b2,c,d,e,dt,t){
         exp(c+d*(t^(exp(a0+a1*dt+a2*dt^2)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))+
         e*(t^(exp(a0+a1*dt+a2*dt^2)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))^2)/
         (1+exp(c+d*(t^(exp(a0+a1*dt+a2*dt^2)-1)*
         (1-t)^(exp(b0+b1*dt+b2*dt^2)-1))+
         e*(t^(exp(a0+a1*dt+a2*dt^2)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))^2))
}

### partial derivative of pi with respect to c_0 (cf., model (2.3))

dpidc0 <- function(ma,c,d,e){
                 exp(c+d*ma+e*ma^2)/(1+exp(c+d*ma+e*ma^2))^2
}

### partial derivative of pi with respect to c_1 (cf., model (2.3))

dpidc1 <- function(ma,c,d,e){
                 ma*exp(c+d*ma+e*ma^2)/(1+exp(c+d*ma+e*ma^2))^2
}

### partial derivative of pi with respect to c_2 (cf., model (2.3))

dpidc2 <- function(ma,c,d,e){
                 ma^2*exp(c+d*ma+e*ma^2)/(1+exp(c+d*ma+e*ma^2))^2
}


####################################################################
###                         Set initial values
####################################################################

beta0 <- c(log(1.5),0,0, log(3),0,0, -3.806156,7.261997,0, 
           -2.094799,4.970579,0, -1.304599,4.698151,0,
           -3.691203,2.896354,0, -4.029930,10.243195,0, 
           -3.252539,4.457725,0, -3.003364,6.902542,0,
           -4.019443,3.803001,0)


####################################################################
###       Set some parameters for the iterative algorithm and 
###       run the algorithm
####################################################################

iter <-0            ### Index of the iterative algorithm
n <-300             ### Maximum number of iterations
itercrit <- 100     ### Criterion for exiting the algorithm

beta1 <- beta0      ### Beta1 denotes the parameter estimates


### The iterative algorithm starts here

while(itercrit > 0.01 & iter < n){

iter <- iter+1
beta0 <- beta1


### Update the five parameters in MA

LM<-matrix(rep(0,25),nc=5)
RM<-rep(0,5)

for (l in 1:sub){

 abcd <- expression(exp(c+d*(t^(exp(a0+a1*dt)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))+
         e*(t^(exp(a0+a1*dt)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))^2)/
         (1+exp(c+d*(t^(exp(a0+a1*dt)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))+
         e*(t^(exp(a0+a1*dt)-1)*(1-t)^(exp(b0+b1*dt+b2*dt^2)-1))^2)))

 (ft1 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[7],d=beta0[8],e=beta0[9], 
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft2 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[10],d=beta0[11],e=beta0[12],
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft3 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[13],d=beta0[14],e=beta0[15],
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft4 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[16],d=beta0[17],e=beta0[18],
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft5 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[19],d=beta0[20],e=beta0[21],
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft6 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[22],d=beta0[23],e=beta0[24],
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft7 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[25],d=beta0[26],e=beta0[27],
  dt=dat[l,9],t=dat[l,10]){} ) )
 (ft8 <- deriv(abcd, c("a0", "a1", "b0", "b1","b2"), 
  function(a0, a1, b0, b1, b2,c=beta0[28],d=beta0[29],e=beta0[30],
  dt=dat[l,9],t=dat[l,10]){} ) )

 qq1<-attr(ft1(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq2<-attr(ft2(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq3<-attr(ft3(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq4<-attr(ft4(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq5<-attr(ft5(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq6<-attr(ft6(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq7<-attr(ft7(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")
 qq8<-attr(ft8(beta0[1],beta0[2],beta0[4],beta0[5],beta0[6]),"gradient")

 qq<-rbind(qq1,qq2,qq3,qq4,qq5,qq6,qq7,qq8)

 pp<-rep(0,8)
 for (i in 1:8){
     pp[i] <- gabcd(beta0[1],beta0[2],beta0[3],beta0[4],beta0[5],beta0[6],
     beta0[i*3+4],beta0[i*3+5],beta0[i*3+6],dat[l,9],dat[l,10])
 }

 LM <- LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
 RM <- RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)

 }

LM <- LM+diag(rep(epsilon,5))

beta1[c(1,2,4,5,6)] <- beta0[c(1,2,4,5,6)]+solve(LM)%*%RM
beta1[3]   <-0


### Update the parameters in the first linkage function

LM<-matrix(rep(0,9),nc=3)
RM<-rep(0,3)

for (l in 1:sub){

 qq1 <- rep(0,3)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[7],beta1[8],beta1[9])
 qq1[2] <- dpidc1(ma,beta1[7],beta1[8],beta1[9])
 qq1[3] <- dpidc2(ma,beta1[7],beta1[8],beta1[9])

 qq<-rbind(qq1,matrix(rep(0,21),nc=3))
 pp<-rep(0,8)

 for (i in 1:8){
     pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
            beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

 LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
 RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)

 }

LM<-LM+diag(rep(epsilon,3))

beta1[7:9]<-beta0[7:9]+solve(LM)%*%RM


### Update the parameters in the second linkage function

LM<-matrix(rep(0,9),nc=3)
RM<-rep(0,3)

for (l in 1:sub){

 qq1 <- rep(0,3)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[10],beta1[11],beta1[12])
 qq1[2] <- dpidc1(ma,beta1[10],beta1[11],beta1[12])
 qq1[3] <- dpidc2(ma,beta1[10],beta1[11],beta1[12])

 qq<-rbind(rep(0,3),qq1,matrix(rep(0,18),nc=3))
 pp<-rep(0,8)

 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

 LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
 RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)

 }

LM<-LM+diag(rep(epsilon,3))

beta1[10:12]<-beta0[10:12]+solve(LM)%*%RM


### Update the parameters in the third linkage function

LM<-matrix(rep(0,9),nc=3)
RM<-rep(0,3)

for (l in 1:sub){

 qq1 <- rep(0,3)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[13],beta1[14],beta1[15])
 qq1[2] <- dpidc1(ma,beta1[13],beta1[14],beta1[15])
 qq1[3] <- dpidc2(ma,beta1[13],beta1[14],beta1[15])

 qq<-rbind(rep(0,3),rep(0,3),qq1,matrix(rep(0,15),nc=3))
 pp<-rep(0,8)
 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)
 }

LM<-LM+diag(rep(epsilon,3))

beta1[13:15]<-beta0[13:15]+solve(LM)%*%RM


### Update the parameters in the fourth linkage function

LM<-matrix(rep(0,4),nc=2)
RM<-rep(0,2)

for (l in 1:sub){

 qq1 <- rep(0,2)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[16],beta1[17],beta1[18])
 qq1[2] <- dpidc1(ma,beta1[16],beta1[17],beta1[18])

 qq<-rbind(matrix(rep(0,6),nc=2),qq1,matrix(rep(0,8),nc=2))

 pp<-rep(0,8)
 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

 LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
 RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)

 }

LM<-LM+diag(rep(epsilon,2))

beta1[16:17]<-beta0[16:17]+solve(LM)%*%RM
beta1[18]   <-0


### Update the parameters in the fifth linkage function

LM<-matrix(rep(0,9),nc=3)
RM<-rep(0,3)

for (l in 1:sub){

 qq1 <- rep(0,3)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[19],beta1[20],beta1[21])
 qq1[2] <- dpidc1(ma,beta1[19],beta1[20],beta1[21])
 qq1[3] <- dpidc2(ma,beta1[19],beta1[20],beta1[21])

 qq<-rbind(matrix(rep(0,12),nc=3),qq1,matrix(rep(0,9),nc=3))

 pp<-rep(0,8)
 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)
 }

LM<-LM+diag(rep(epsilon,3))
beta1[19:21]<-beta0[19:21]+solve(LM)%*%RM


### Update the parameters in the sixth linkage function

LM<-matrix(rep(0,4),nc=2)
RM<-rep(0,2)

for (l in 1:sub){

 qq1 <- rep(0,2)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[22],beta1[23],beta1[24])
 qq1[2] <- dpidc1(ma,beta1[22],beta1[23],beta1[24])

 qq<-rbind(matrix(rep(0,10),nc=2),qq1,matrix(rep(0,4),nc=2))

 pp<-rep(0,8)
 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

 LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
 RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)
 }

LM<-LM+diag(rep(epsilon,2))

beta1[22:23]<-beta0[22:23]+solve(LM)%*%RM
beta1[24]   <-0


### Update the parameters in the seventh linkage function

LM<-matrix(rep(0,9),nc=3)
RM<-rep(0,3)

for (l in 1:sub){

 qq1 <- rep(0,3)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[25],beta1[26],beta1[27])
 qq1[2] <- dpidc1(ma,beta1[25],beta1[26],beta1[27])
 qq1[3] <- dpidc2(ma,beta1[25],beta1[26],beta1[27])

 qq<-rbind(matrix(rep(0,18),nc=3),qq1,rep(0,3))

 pp<-rep(0,8)
 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)
 }

LM<-LM+diag(rep(epsilon,3))
beta1[25:27]<-beta0[25:27]+solve(LM)%*%RM


### Update the parameters in the eighth linkage function

LM<-matrix(rep(0,9),nc=3)
RM<-rep(0,3)

for (l in 1:sub){

 qq1 <- rep(0,3)

 dt=dat[l,9]
 t=dat[l,10]

 a=beta1[1]+beta1[2]*dt+beta1[3]*dt^2
 b=beta1[4]+beta1[5]*dt+beta1[6]*dt^2

 ma  <- t^(exp(a)-1)*(1-t)^(exp(b)-1)

 qq1[1] <- dpidc0(ma,beta1[28],beta1[29],beta1[30])
 qq1[2] <- dpidc1(ma,beta1[28],beta1[29],beta1[30])
 qq1[3] <- dpidc2(ma,beta1[28],beta1[29],beta1[30])

 qq<-rbind(matrix(rep(0,21),nc=3),qq1)

 pp<-rep(0,8)
 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

 LM<-LM+t(qq)%*%diag(1/(pp*(1-pp)))%*%qq
 RM<-RM+t(qq)%*%diag(1/(pp*(1-pp)))%*%(dat[l,1:8]-pp)

 }

LM<-LM+diag(rep(epsilon,3))
beta1[28:30]<-beta0[28:30]+solve(LM)%*%RM

### Update the quiting criterion value

itercrit <- sqrt(sum((beta1-beta0)^2))/sqrt(sum(beta1^2))

### Write results to outfile

write(c(iter,itercrit),outfile)
write(beta1,outfile)
write('',outfile)

### Compute the QIC criterion value

qic <- 0

for (l in 1:sub){

 pp<-rep(0,8)

 for (i in 1:8){
 pp[i]<-gabcd(beta1[1],beta1[2],beta1[3],beta1[4],beta1[5],beta1[6],
        beta1[i*3+4],beta1[i*3+5],beta1[i*3+6],dat[l,9],dat[l,10])
 }

 qic <- qic + sum(dat[l,1:8]*log(pp/(1-pp))+log(1-pp))

 }

qic <- -2*qic+2*27

write(qic,outfile)

}