# -----------------------------------------------------------
# RPD analysis
# -----------------------------------------------------------

# -----------------------------------------------------------
# read data
# SET WORKING DIRECTORY TO READ EXAMPLE DATA
setwd("c:/WHERE IS YOUR DATA/data")
qscores.dat <- read.table("quality_data.txt",
      col.names=c("Pot","Plant","Year","Time","Week","Q","Temp","Light","Water","Cultivar",
                  "Nitrogen","Potassium","Chain","Q10","Q4"),header=TRUE)
# -----------------------------------------------------------


# -----------------------------------------------------------
# set factor levels
qscores.dat$Temp <- factor(qscores.dat$Temp,labels=c("16C","21C"))
qscores.dat$Light <- factor(qscores.dat$Light,labels=c("300lux","600lux"))
qscores.dat$Water <- factor(qscores.dat$Water,labels=c("Control","Fluct"))
qscores.dat$Cultivar <- factor(qscores.dat$Cultivar,labels=c("Batik","Mariette"))
qscores.dat$Nitrogen <- factor(qscores.dat$Nitrogen,labels=c("Low","High"))
qscores.dat$Potassium <- factor(qscores.dat$Potassium,labels=c("Low","High"))
qscores.dat$Chain <- factor(qscores.dat$Chain,labels=c("Mini","Supm"))
# -----------------------------------------------------------


# -----------------------------------------------------------
# source repolr and gee library
library(gee)
# SET SOURCE FILE TO FIND FUNCTION FOR REPOLR
source("h:/Nick/papers/JAS/analysis/repolr/repolr_functions.R")
# -----------------------------------------------------------



# -----------------------------------------------------------
# fit overall model
overall_mod <- repolr(Q4~Time,qscores.dat,subjects="Plant",categories=4,
                        times=c(1,2,3,4,5,6,7),corstr="ar1")
# -----------------------------------------------------------


# -----------------------------------------------------------
# fit model for each plant i = 1...128
toler <- c(1,0.75,0.5,0.05,0.001,0.0001)
planti <- vector(length=128)
oplanti <- planti
for (i in 1:128){
 qscoresi <- data.frame(qscores.dat$Q4,qscores.dat$Plant,
          qscores.dat$Time,as.numeric(qscores.dat$Plant!=i))
 colnames(qscoresi) <- c("Q4","Plant","Week","Pick")
 datai <- ord.expand(scores="Q4",data=qscoresi,subjects="Plant",categories=4)$exdata
 datai$cuts <- as.numeric(datai$cuts)
 datai$cuts[datai$cuts==1] <- overall_mod[["gee"]]$coeff[1]
 datai$cuts[datai$cuts==2] <- overall_mod[["gee"]]$coeff[2]
 datai$cuts[datai$cuts==3] <- overall_mod[["gee"]]$coeff[3]
 picker <- gee(Q4~factor(Pick):Week-1+offset(cuts),datai,id=subjects,
        family="binomial",corstr="fixed",scale.value=1,scale.fix=TRUE,
        b=c(overall_mod[["gee"]]$coeff[4],overall_mod[["gee"]]$coeff[4]),
        R=overall_mod[["gee"]]$working.correlation,tol=toler[3])
 planti[i] <- picker$coeff[1]
 oplanti[i] <- picker$coeff[2]
}
# -----------------------------------------------------------

# -----------------------------------------------------------
# prediction function 1
pred_funct1 <- function(coeff,week){
 pred1 <- exp(overall_mod[["gee"]]$coeff[1]+coeff*week)/(1+exp(overall_mod[["gee"]]$coeff[1]+coeff*week))
 pred2 <- exp(overall_mod[["gee"]]$coeff[2]+coeff*week)/(1+exp(overall_mod[["gee"]]$coeff[2]+coeff*week))
 pred3 <- exp(overall_mod[["gee"]]$coeff[3]+coeff*week)/(1+exp(overall_mod[["gee"]]$coeff[3]+coeff*week))
 prob1 <- pred1
 prob2 <- pred2-pred1
 prob3 <- pred3-pred2
 prob4 <- 1-pred3
 Omean <- 4*prob4+3*prob3+2*prob2+1*prob1
 list(Omean=Omean,prob1=prob1,prob2=prob2,prob3=prob3,prob4=prob4)
}
# -----------------------------------------------------------

# -----------------------------------------------------------
# prediction function 2
pred_funct2 <- function(coeff,week){
 prob <- exp(-coeff*week)/(1+exp(-coeff*week))
 list(prob=prob)
}
# -----------------------------------------------------------

# -----------------------------------------------------------
# inverse prediction function 2
ipred_funct2 <- function(prob,coeff){
 week <- 3-log(prob/(1-prob))/coeff
 list(week=week)
}
# -----------------------------------------------------------

# -----------------------------------------------------------
# Fitted model summary
Rqscores.dat <- data.frame(qscores.dat$Plant[seq(1,896,7)],qscores.dat$Temp[seq(1,896,7)],
                           qscores.dat$Light[seq(1,896,7)],qscores.dat$Water[seq(1,896,7)],
                           qscores.dat$Cultivar[seq(1,896,7)],qscores.dat$Nitrogen[seq(1,896,7)],
                           qscores.dat$Potassium[seq(1,896,7)],qscores.dat$Chain[seq(1,896,7)],
                           planti,oplanti)
colnames(Rqscores.dat) <- c("Plant","Temp","Light","Water","Cultivar","Nitrogen",
                              "Potassium","Chain","Coeff1","Coeff2")
# means
mean.tab <- tapply(Rqscores.dat$Coeff1,list(Rqscores.dat$Chain,Rqscores.dat$Cultivar,
                          Rqscores.dat$Potassium,Rqscores.dat$Nitrogen),mean,na.rm=TRUE)
qsvkn_1111 <- mean.tab["Mini","Batik","Low","Low"]
qsvkn_1112 <- mean.tab["Mini","Batik","Low","High"]
qsvkn_1121 <- mean.tab["Mini","Batik","High","Low"]
qsvkn_1122 <- mean.tab["Mini","Batik","High","High"]
qsvkn_1211 <- mean.tab["Mini","Mariette","Low","Low"]
qsvkn_1212 <- mean.tab["Mini","Mariette","Low","High"]
qsvkn_1221 <- mean.tab["Mini","Mariette","High","Low"]
qsvkn_1222 <- mean.tab["Mini","Mariette","High","High"]
qsvkn_2111 <- mean.tab["Supm","Batik","Low","Low"]
qsvkn_2112 <- mean.tab["Supm","Batik","Low","High"]
qsvkn_2121 <- mean.tab["Supm","Batik","High","Low"]
qsvkn_2122 <- mean.tab["Supm","Batik","High","High"]
qsvkn_2211 <- mean.tab["Supm","Mariette","Low","Low"]
qsvkn_2212 <- mean.tab["Supm","Mariette","Low","High"]
qsvkn_2221 <- mean.tab["Supm","Mariette","High","Low"]
qsvkn_2222 <- mean.tab["Supm","Mariette","High","High"]
Rqscores.dat
# maximum
max.tab <- tapply(Rqscores.dat$Coeff1,list(Rqscores.dat$Chain,Rqscores.dat$Cultivar,
                          Rqscores.dat$Potassium,Rqscores.dat$Nitrogen),max,na.rm=TRUE)
usvkn_1111 <- max.tab["Mini","Batik","Low","Low"]
usvkn_1112 <- max.tab["Mini","Batik","Low","High"]
usvkn_1121 <- max.tab["Mini","Batik","High","Low"]
usvkn_1122 <- max.tab["Mini","Batik","High","High"]
usvkn_1211 <- max.tab["Mini","Mariette","Low","Low"]
usvkn_1212 <- max.tab["Mini","Mariette","Low","High"]
usvkn_1221 <- max.tab["Mini","Mariette","High","Low"]
usvkn_1222 <- max.tab["Mini","Mariette","High","High"]
usvkn_2111 <- max.tab["Supm","Batik","Low","Low"]
usvkn_2112 <- max.tab["Supm","Batik","Low","High"]
usvkn_2121 <- max.tab["Supm","Batik","High","Low"]
usvkn_2122 <- max.tab["Supm","Batik","High","High"]
usvkn_2211 <- max.tab["Supm","Mariette","Low","Low"]
usvkn_2212 <- max.tab["Supm","Mariette","Low","High"]
usvkn_2221 <- max.tab["Supm","Mariette","High","Low"]
usvkn_2222 <- max.tab["Supm","Mariette","High","High"]
# lower
min.tab <- tapply(Rqscores.dat$Coeff1,list(Rqscores.dat$Chain,Rqscores.dat$Cultivar,
                          Rqscores.dat$Potassium,Rqscores.dat$Nitrogen),min,na.rm=TRUE)
lsvkn_1111 <- min.tab["Mini","Batik","Low","Low"]
lsvkn_1112 <- min.tab["Mini","Batik","Low","High"]
lsvkn_1121 <- min.tab["Mini","Batik","High","Low"]
lsvkn_1122 <- min.tab["Mini","Batik","High","High"]
lsvkn_1211 <- min.tab["Mini","Mariette","Low","Low"]
lsvkn_1212 <- min.tab["Mini","Mariette","Low","High"]
lsvkn_1221 <- min.tab["Mini","Mariette","High","Low"]
lsvkn_1222 <- min.tab["Mini","Mariette","High","High"]
lsvkn_2111 <- min.tab["Supm","Batik","Low","Low"]
lsvkn_2112 <- min.tab["Supm","Batik","Low","High"]
lsvkn_2121 <- min.tab["Supm","Batik","High","Low"]
lsvkn_2122 <- min.tab["Supm","Batik","High","High"]
lsvkn_2211 <- min.tab["Supm","Mariette","Low","Low"]
lsvkn_2212 <- min.tab["Supm","Mariette","Low","High"]
lsvkn_2221 <- min.tab["Supm","Mariette","High","Low"]
lsvkn_2222 <- min.tab["Supm","Mariette","High","High"]
# -----------------------------------------------------------


# -----------------------------------------------------------
# Kappa agreement statistic
# -----------------------------------------------------------
predcat <- NULL
for (i in 1:128){
 expect <- pred_funct1(planti[i],c(-3,-2,-1,0,1,2,3))
 probs <- matrix(c(expect$prob1,expect$prob2,expect$prob3,expect$prob4),
                           nrow=4,ncol=7,byrow=TRUE)
 expcat <- vector(length=7)
 for (j in 1:7){
  expcat[j] <- which.max(probs[,j])
 }
 predcat <- append(predcat,expcat)
}
library(psy)
nqscores.dat <- data.frame(qscores.dat,predcat)
Ksvkn_1111 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_1112 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="High",]
Ksvkn_1121 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_1122 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="High",]
Ksvkn_1211 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_1212 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="High",]
Ksvkn_1221 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_1222 <- nqscores.dat[nqscores.dat$Chain=="Mini" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="High",]
Ksvkn_2111 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_2112 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="High",]
Ksvkn_2121 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_2122 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Batik" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="High",]
Ksvkn_2211 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_2212 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="Low" & nqscores.dat$Nitrogen=="High",]
Ksvkn_2221 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="Low",]
Ksvkn_2222 <- nqscores.dat[nqscores.dat$Chain=="Supm" & nqscores.dat$Cultivar=="Mariette" &
                           nqscores.dat$Potassium=="High" & nqscores.dat$Nitrogen=="High",]
amean_data <- list(Ksvkn_1111,Ksvkn_1112,Ksvkn_1121,Ksvkn_1122,
                   Ksvkn_1211,Ksvkn_1212,Ksvkn_1221,Ksvkn_1222,
                   Ksvkn_2111,Ksvkn_2112,Ksvkn_2121,Ksvkn_2122,
                   Ksvkn_2211,Ksvkn_2212,Ksvkn_2221,Ksvkn_2222)
svkn_kappa <- vector(length=16)
for (k in 1:16){
 svkn_kappa[k] <- wkappa(matrix(c(amean_data[[k]]$Q4,amean_data[[k]]$predcat),
                   ncol=2,nrow=56,byrow=FALSE),weights="squared")$kappa
}
# -----------------------------------------------------------


# -----------------------------------------------------------
# Observed data summary
count1 <- function(x){
 length(x[x==1])
}
count2 <- function(x){
 length(x[x==2])
}
count3 <- function(x){
 length(x[x==3])
}
count4 <- function(x){
 length(x[x==4])
}
count1.tab <- tapply(qscores.dat$Q4,list(qscores.dat$Chain,qscores.dat$Cultivar,
                          qscores.dat$Potassium,qscores.dat$Nitrogen,
                          qscores.dat$Week),count1)
count2.tab <- tapply(qscores.dat$Q4,list(qscores.dat$Chain,qscores.dat$Cultivar,
                          qscores.dat$Potassium,qscores.dat$Nitrogen,
                          qscores.dat$Week),count2)
count3.tab <- tapply(qscores.dat$Q4,list(qscores.dat$Chain,qscores.dat$Cultivar,
                          qscores.dat$Potassium,qscores.dat$Nitrogen,
                          qscores.dat$Week),count3)
count4.tab <- tapply(qscores.dat$Q4,list(qscores.dat$Chain,qscores.dat$Cultivar,
                          qscores.dat$Potassium,qscores.dat$Nitrogen,
                          qscores.dat$Week),count4)
# count 1
m1svkn_1111 <- count1.tab["Mini","Batik","Low","Low",]
m1svkn_1112 <- count1.tab["Mini","Batik","Low","High",]
m1svkn_1121 <- count1.tab["Mini","Batik","High","Low",]
m1svkn_1122 <- count1.tab["Mini","Batik","High","High",]
m1svkn_1211 <- count1.tab["Mini","Mariette","Low","Low",]
m1svkn_1212 <- count1.tab["Mini","Mariette","Low","High",]
m1svkn_1221 <- count1.tab["Mini","Mariette","High","Low",]
m1svkn_1222 <- count1.tab["Mini","Mariette","High","High",]
m1svkn_2111 <- count1.tab["Supm","Batik","Low","Low",]
m1svkn_2112 <- count1.tab["Supm","Batik","Low","High",]
m1svkn_2121 <- count1.tab["Supm","Batik","High","Low",]
m1svkn_2122 <- count1.tab["Supm","Batik","High","High",]
m1svkn_2211 <- count1.tab["Supm","Mariette","Low","Low",]
m1svkn_2212 <- count1.tab["Supm","Mariette","Low","High",]
m1svkn_2221 <- count1.tab["Supm","Mariette","High","Low",]
m1svkn_2222 <- count1.tab["Supm","Mariette","High","High",]
# count 2
m2svkn_1111 <- count2.tab["Mini","Batik","Low","Low",]
m2svkn_1112 <- count2.tab["Mini","Batik","Low","High",]
m2svkn_1121 <- count2.tab["Mini","Batik","High","Low",]
m2svkn_1122 <- count2.tab["Mini","Batik","High","High",]
m2svkn_1211 <- count2.tab["Mini","Mariette","Low","Low",]
m2svkn_1212 <- count2.tab["Mini","Mariette","Low","High",]
m2svkn_1221 <- count2.tab["Mini","Mariette","High","Low",]
m2svkn_1222 <- count2.tab["Mini","Mariette","High","High",]
m2svkn_2111 <- count2.tab["Supm","Batik","Low","Low",]
m2svkn_2112 <- count2.tab["Supm","Batik","Low","High",]
m2svkn_2121 <- count2.tab["Supm","Batik","High","Low",]
m2svkn_2122 <- count2.tab["Supm","Batik","High","High",]
m2svkn_2211 <- count2.tab["Supm","Mariette","Low","Low",]
m2svkn_2212 <- count2.tab["Supm","Mariette","Low","High",]
m2svkn_2221 <- count2.tab["Supm","Mariette","High","Low",]
m2svkn_2222 <- count2.tab["Supm","Mariette","High","High",]
# count 3
m3svkn_1111 <- count3.tab["Mini","Batik","Low","Low",]
m3svkn_1112 <- count3.tab["Mini","Batik","Low","High",]
m3svkn_1121 <- count3.tab["Mini","Batik","High","Low",]
m3svkn_1122 <- count3.tab["Mini","Batik","High","High",]
m3svkn_1211 <- count3.tab["Mini","Mariette","Low","Low",]
m3svkn_1212 <- count3.tab["Mini","Mariette","Low","High",]
m3svkn_1221 <- count3.tab["Mini","Mariette","High","Low",]
m3svkn_1222 <- count3.tab["Mini","Mariette","High","High",]
m3svkn_2111 <- count3.tab["Supm","Batik","Low","Low",]
m3svkn_2112 <- count3.tab["Supm","Batik","Low","High",]
m3svkn_2121 <- count3.tab["Supm","Batik","High","Low",]
m3svkn_2122 <- count3.tab["Supm","Batik","High","High",]
m3svkn_2211 <- count3.tab["Supm","Mariette","Low","Low",]
m3svkn_2212 <- count3.tab["Supm","Mariette","Low","High",]
m3svkn_2221 <- count3.tab["Supm","Mariette","High","Low",]
m3svkn_2222 <- count3.tab["Supm","Mariette","High","High",]
# count 4
m4svkn_1111 <- count4.tab["Mini","Batik","Low","Low",]
m4svkn_1112 <- count4.tab["Mini","Batik","Low","High",]
m4svkn_1121 <- count4.tab["Mini","Batik","High","Low",]
m4svkn_1122 <- count4.tab["Mini","Batik","High","High",]
m4svkn_1211 <- count4.tab["Mini","Mariette","Low","Low",]
m4svkn_1212 <- count4.tab["Mini","Mariette","Low","High",]
m4svkn_1221 <- count4.tab["Mini","Mariette","High","Low",]
m4svkn_1222 <- count4.tab["Mini","Mariette","High","High",]
m4svkn_2111 <- count4.tab["Supm","Batik","Low","Low",]
m4svkn_2112 <- count4.tab["Supm","Batik","Low","High",]
m4svkn_2121 <- count4.tab["Supm","Batik","High","Low",]
m4svkn_2122 <- count4.tab["Supm","Batik","High","High",]
m4svkn_2211 <- count4.tab["Supm","Mariette","Low","Low",]
m4svkn_2212 <- count4.tab["Supm","Mariette","Low","High",]
m4svkn_2221 <- count4.tab["Supm","Mariette","High","Low",]
m4svkn_2222 <- count4.tab["Supm","Mariette","High","High",]
# -----------------------------------------------------------


# -----------------------------------------------------------
# all data
amean_data <- list(qsvkn_1111,qsvkn_1112,qsvkn_1121,qsvkn_1122,
                   qsvkn_1211,qsvkn_1212,qsvkn_1221,qsvkn_1222,
                   qsvkn_2111,qsvkn_2112,qsvkn_2121,qsvkn_2122,
                   qsvkn_2211,qsvkn_2212,qsvkn_2221,qsvkn_2222)
amax_data <- list(usvkn_1111,usvkn_1112,usvkn_1121,usvkn_1122,
                   usvkn_1211,usvkn_1212,usvkn_1221,usvkn_1222,
                   usvkn_2111,usvkn_2112,usvkn_2121,usvkn_2122,
                   usvkn_2211,usvkn_2212,usvkn_2221,usvkn_2222)
amin_data <- list(lsvkn_1111,lsvkn_1112,lsvkn_1121,lsvkn_1122,
                   lsvkn_1211,lsvkn_1212,lsvkn_1221,lsvkn_1222,
                   lsvkn_2111,lsvkn_2112,lsvkn_2121,lsvkn_2122,
                   lsvkn_2211,lsvkn_2212,lsvkn_2221,lsvkn_2222)
count1_data <- list(m1svkn_1111,m1svkn_1112,m1svkn_1121,m1svkn_1122,
                   m1svkn_1211,m1svkn_1212,m1svkn_1221,m1svkn_1222,
                   m1svkn_2111,m1svkn_2112,m1svkn_2121,m1svkn_2122,
                   m1svkn_2211,m1svkn_2212,m1svkn_2221,m1svkn_2222)
count2_data <- list(m2svkn_1111,m2svkn_1112,m2svkn_1121,m2svkn_1122,
                   m2svkn_1211,m2svkn_1212,m2svkn_1221,m2svkn_1222,
                   m2svkn_2111,m2svkn_2112,m2svkn_2121,m2svkn_2122,
                   m2svkn_2211,m2svkn_2212,m2svkn_2221,m2svkn_2222)
count3_data <- list(m3svkn_1111,m3svkn_1112,m3svkn_1121,m3svkn_1122,
                   m3svkn_1211,m3svkn_1212,m3svkn_1221,m3svkn_1222,
                   m3svkn_2111,m3svkn_2112,m3svkn_2121,m3svkn_2122,
                   m3svkn_2211,m3svkn_2212,m3svkn_2221,m3svkn_2222)
count4_data <- list(m4svkn_1111,m4svkn_1112,m4svkn_1121,m4svkn_1122,
                   m4svkn_1211,m4svkn_1212,m4svkn_1221,m4svkn_1222,
                   m4svkn_2111,m4svkn_2112,m4svkn_2121,m4svkn_2122,
                   m4svkn_2211,m4svkn_2212,m4svkn_2221,m4svkn_2222)
# all labels
alab_data <- list(c(1,1,1,1),c(1,1,1,2),c(1,1,2,1),c(1,1,2,2),
                  c(1,2,1,1),c(1,2,1,2),c(1,2,2,1),c(1,2,2,2),
                  c(2,1,1,1),c(2,1,1,2),c(2,1,2,1),c(2,1,2,2),
                  c(2,2,1,1),c(2,2,1,2),c(2,2,2,1),c(2,2,2,2))
# -----------------------------------------------------------



# -----------------------------------------------------------
# Figure 1
# -----------------------------------------------------------

# -----------------------------------------------------------
# plot layout
X11()
nf <- layout(matrix(c(1,2,3,4,5,
                      1,6,7,8,9,
                      1,10,11,12,13,
                      1,14,15,16,17,
                      0,18,18,18,18),5,5,byrow=TRUE),
                      c(1,4,4,4,4), c(4,4,4,4,1), TRUE)
layout.show(nf)
# -----------------------------------------------------------

# -----------------------------------------------------------
# y label
par(font.axis=1,font.lab=1,font.main=1,pty="m",mai=c(0,0,0,0),
       bty="n",xaxt="n",yaxt="n")
plot(x=1:2,y=1:2,pch=" ",ylim=c(1,2),xlim=c(1,2))
text(x=1.3,y=1.5,"Quality Score",srt=90,cex=2)
text(x=c(rep(1.85,4)),y=c(1.218,1.187,1.010,0.975),
                 c(4,3,2,1),cex=1.6)
text(x=c(rep(1.85,4)),y=c(1.218,1.187,1.010,0.975)+0.27,
                 c(4,3,2,1),cex=1.6)
text(x=c(rep(1.85,4)),y=c(1.218,1.187,1.010,0.975)+0.54,
                 c(4,3,2,1),cex=1.6)
text(x=c(rep(1.85,4)),y=c(1.218,1.187,1.010,0.975)+0.81,
                 c(4,3,2,1),cex=1.6)
# -----------------------------------------------------------

# -----------------------------------------------------------
# functions
xtcks_fun <- function(x){
 for(i in 1:x){
  j <- i-1
  lines(x=c(j,j),y=c(-0.05,0.00),lty=1,lwd=1)
 }
}

yticks_fun <- function(x,ytick){
 for(i in 1:x){
  lines(x=c(-0.9,-0.5),y=c(ytick[i],ytick[i]),lty=1,lwd=1)
 }
}
error_bars <- function(lower,upper){
 for (k in 1:7){
  j <- k-1
  arrows(x0=j,x1=j,y0=lower[k],y1=upper[k],
         angle=90,code=3,length=0.03,lty=2)
 }
}
lab_fun <- function(lab){
 text(x=c(4.8,5.2,5.6,6.0),y=c(1.00,1.00,1.00,1.00),
                                 c("S","V","K","N"),cex=1.1)
 text(x=c(4.8,5.2,5.6,6.0),y=c(0.93,0.93,0.93,0.93),
                                 lab,cex=1.1)
}
# -----------------------------------------------------------



# -----------------------------------------------------------
# plots
# point sizes for plotting
p_size <- c(0,0.5,0.7,0.9,1.1,1.3,1.5,1.7,1.9)
# location of scores 4, 3, 2 and 1 on interval [0,1]
# determined numerically
ord_y <- c(0.00,0.15,0.86,1.00)
for (p in 1:16){
 par(font.axis=1,font.lab=1,font.main=1,pty="s",mai=c(0,0,0,0),
       bty="o",xaxt="n",yaxt="n",tcl=0.3)
 plot(x=seq(0,6,by=0.05),y=pred_funct2(amean_data[[p]],seq(-3,3,by=0.05))$prob,lwd=2,
         las=1,cex.axis=1.3,cex.lab=1.3,cex=1.3,xlab="Time (Years)",ylim=c(0,1),
             xlim=c(-0.5,6.5),type="l")
 lines(x=seq(0,6,by=0.05),y=pred_funct2(amin_data[[p]],seq(-3,3,by=0.05))$prob,lwd=1,lty=2)
 lines(x=seq(0,6,by=0.05),y=pred_funct2(amax_data[[p]],seq(-3,3,by=0.05))$prob,lwd=1,lty=2)
 xtcks_fun(7)
 yticks_fun(x=4,ytick=ord_y)
 lab_fun(alab_data[[p]])
 for (k in 1:7){
  sym_size1 <- p_size[count1_data[[p]][k]+1]
  sym_size2 <- p_size[count2_data[[p]][k]+1]
  sym_size3 <- p_size[count3_data[[p]][k]+1]
  sym_size4 <- p_size[count4_data[[p]][k]+1]
  points(x=(k-1),y=ord_y[4],cex=sym_size4,pch=19)
  points(x=(k-1),y=ord_y[3],cex=sym_size3,pch=19)
  points(x=(k-1),y=ord_y[2],cex=sym_size2,pch=19)
  points(x=(k-1),y=ord_y[1],cex=sym_size1,pch=19)
 }
 kap_tit <- expression(italic(kappa)==phantom(1))
 text(x=0.5,y=0.095,kap_tit,cex=1.9)
 text(x=1.5,y=0.085,as.character(format(round(svkn_kappa[p],2),digits=3,trim=TRUE,
       justify="right",scientific=FALSE,nsmall=2,width=3)),cex=1.2)
}
# -----------------------------------------------------------


# -----------------------------------------------------------
# x label
par(font.axis=1,font.lab=1,font.main=1,pty="m",mai=c(0,0,0,0),
       bty="n",xaxt="n",yaxt="n")
plot(x=1:2,y=1:2,pch=" ",ylim=c(1,2),xlim=c(1,2))
text(x=1.5,y=1.3,"Time in Home-life (weeks)",srt=0,cex=2)
text(x=c(0.99,1.025,1.061,1.097,1.133,1.169,1.205),
                y=rep(1.82,7),c(0,1,2,3,4,5,6),cex=1.6)
text(x=c(0.99,1.025,1.061,1.097,1.133,1.169,1.205)+0.27,
                y=rep(1.82,7),c(0,1,2,3,4,5,6),cex=1.6)
text(x=c(0.99,1.025,1.061,1.097,1.133,1.169,1.205)+0.54,
                y=rep(1.82,7),c(0,1,2,3,4,5,6),cex=1.6)
text(x=c(0.99,1.025,1.061,1.097,1.133,1.169,1.205)+0.81,
                y=rep(1.82,7),c(0,1,2,3,4,5,6),cex=1.6)
# -----------------------------------------------------------



# -----------------------------------------------------------
# Modelling
# -----------------------------------------------------------


# -----------------------------------------------------------
# set-up factors and contrasts
# -----------------------------------------------------------

# factors
TT <- Rqscores.dat$Temp; RR <- Rqscores.dat$Light; WW <- Rqscores.dat$Water
SS <- Rqscores.dat$Chain; VV <- Rqscores.dat$Cultivar; NN <- Rqscores.dat$Nitrogen
KK <- Rqscores.dat$Potassium; PP <- Rqscores.dat$Plant

# contrasts
contrasts(TT) <- contr.treatment(2,base=1); attributes(contrasts(TT))$dimnames[[2]] <- "T+"
contrasts(RR) <- contr.treatment(2,base=1); attributes(contrasts(RR))$dimnames[[2]] <- "R+"
contrasts(WW) <- contr.treatment(2,base=1); attributes(contrasts(WW))$dimnames[[2]] <- "W+"
contrasts(SS) <- contr.treatment(2,base=1); attributes(contrasts(SS))$dimnames[[2]] <- "S+"
contrasts(VV) <- contr.treatment(2,base=1); attributes(contrasts(VV))$dimnames[[2]] <- "V+"
contrasts(NN) <- contr.treatment(2,base=1); attributes(contrasts(NN))$dimnames[[2]] <- "N+"
contrasts(KK) <- contr.treatment(2,base=1); attributes(contrasts(KK))$dimnames[[2]] <- "K+"


#------------------------------------------------------------
# Joint model analysis
#------------------------------------------------------------

# response variate
rate <- planti

# ANOVA for rate
rate_aov <- aov(rate~NN+SS+VV+TT+RR+KK+WW+TT:RR+TT:VV+RR:VV+TT:RR:VV+VV:TT:RR)

# initial settings
mwt <- c(rep(1,128)); nunits <- length(rate)
ctiter <- 0; ogstop <- 1
mdeviance <- NULL; mresiddf <- NULL; ddeviance <- NULL; dresiddf <- NULL
mddf <- NULL; EQD <- NULL; EQU <- NULL; APL <- NULL

# while loop
while (ogstop>0.01 & ctiter<=5){

 # start model fitting

 # mean model
 mean_mod <- glm(rate~NN+VV+SS+TT+RR+KK+WW+TT:RR+TT:VV+RR:VV+VV:TT:RR+RR:KK+WW:TT,
             family=gaussian(link="identity"),weight=mwt, x=TRUE)
 sum_mean <- summary.glm(mean_mod)
 mlev <- hatvalues(mean_mod)
 mfit <- fitted.values(mean_mod)
 mdeviance <- append(mdeviance,sum_mean$deviance)
 mresiddf <- append(mresiddf,sum_mean$df.residual)

 # modify effect of large leverages
 dispwt <- (mlev>1)*0.5+(mlev<=1)*mlev
 dispwt <- 1-dispwt

 # response for dispersion model
 disp_resp <- (mfit-rate)**2
 disp_resp <- disp_resp/dispwt

 # dispersion model
 disp_mod <- glm(disp_resp~WW, family=Gamma(link="log"),weight=dispwt)
 sum_disp <- summary.glm(disp_mod, dispersion=2)
 dfit <- fitted.values(disp_mod)
 ddeviance <- append(ddeviance,sum_disp$deviance)
 dresiddf <- append(dresiddf,sum_disp$df.residual)

 # update weights for mean model
 mwt <- 1/dfit

 # statistics 
 cphi <- dfit
 mddf <- append(mddf,nunits-(nunits-sum_mean$df.residual)-(nunits-sum_disp$df.residual))
 
 # EQD: extended quasi-deviance (with h-adjustment)
 mcdev <- ((mfit-rate)**2)/(1-mlev)
 cEQD <- (mcdev/cphi)+log(2*pi*cphi)
 EQD <- append(EQD,sum(cEQD))

 # EQU: extended quasi-deviance (without h-adjustment)
 mdev <- ((mfit-rate)**2)
 cEQU <- (mdev/cphi)+log(2*pi*cphi)
 EQU <- append(EQU,sum(cEQU))

 # APL: adjusted profile likelihood
 wstar <- 1/dfit
 imm <- solve(t(mean_mod$x)%*%diag(wstar)%*%mean_mod$x)
 APL <- append(APL,sum(cEQU)-sum(log(eigen(imm)$values*2*pi)))

 # monitor convergence
 ctiter <- ctiter+1
 newcoeff <- c(mean_mod$coeff,disp_mod$coeff)
 if (ctiter>1) {crit <- 100*abs(1-sqrt(sum((fullcoeff/newcoeff)^2)/length(fullcoeff)))}
 fullcoeff <- newcoeff
 if (ctiter<=2){ogstop <- 1} else {ogstop <- crit}
 
 # monitor data
 md_monitor <- matrix(0,nrow=ctiter,ncol=8,dimnames = list(1:ctiter,
       c("m_dev","m_df", "d_dev","d_df","EQD","EQU","APL","df")))
 md_monitor[,1] <- mdeviance
 md_monitor[,2] <- mresiddf
 md_monitor[,3] <- ddeviance
 md_monitor[,4] <- dresiddf
 md_monitor[,5] <- EQD
 md_monitor[,6] <- EQU
 md_monitor[,7] <- APL
 md_monitor[,8] <- mddf

 # output
 list(monitor=md_monitor,coeff=fullcoeff)

} # end while


#------------------------------------------------------------
# End of joint model analysis
#-----------------------------------------------------------


#------------------------------------------------------------
# Model summary
#------------------------------------------------------------
sum_disp
sum_mean
md_monitor
# -----------------------------------------------------------


#------------------------------------------------------------
# Summary Plots
#------------------------------------------------------------

# -----------------------------------------------------------
# Figure 3
# -----------------------------------------------------------

# residuals Normal order statistics
X11()
par(font.axis=1,font.lab=1,font.main=1,pty="s")
qqnorm(residuals(mean_mod,type=c("response")),pch="+",cex.axis=1.3,cex.lab=1.3,cex=1.3,
         xlab="Normal Order Statistics",ylim=c(-0.8,0.8),xlim=c(-3,3),
         main="",lab=c(5,8,5),las=1,ylab="Residuals")
qqline(residuals(mean_mod,type=c("response")),lty=2)
lines(-4:4,rep(0,9),type="l")
lines(rep(0,9),-4:4,type="l")

# -----------------------------------------------------------
# Figure 2
# -----------------------------------------------------------

# residuals to show importance of dispersion model
disp_resp <- rate-mfit
disp_resp2 <- disp_resp^2
svkn_ind <- vector(length=128)
PPsvkn_1111 <- PP[SS=="Mini" & VV=="Batik" & KK=="Low" & NN=="Low"]
PPsvkn_1112 <- PP[SS=="Mini" & VV=="Batik" & KK=="Low" & NN=="High"]
PPsvkn_1121 <- PP[SS=="Mini" & VV=="Batik" & KK=="High" & NN=="Low"]
PPsvkn_1122 <- PP[SS=="Mini" & VV=="Batik" & KK=="High" & NN=="High"]
PPsvkn_1211 <- PP[SS=="Mini" & VV=="Mariette" & KK=="Low" & NN=="Low"]
PPsvkn_1212 <- PP[SS=="Mini" & VV=="Mariette" & KK=="Low" & NN=="High"]
PPsvkn_1221 <- PP[SS=="Mini" & VV=="Mariette" & KK=="High" & NN=="Low"]
PPsvkn_1222 <- PP[SS=="Mini" & VV=="Mariette" & KK=="High" & NN=="High"]
PPsvkn_2111 <- PP[SS=="Supm" & VV=="Batik" & KK=="Low" & NN=="Low"]
PPsvkn_2112 <- PP[SS=="Supm" & VV=="Batik" & KK=="Low" & NN=="High"]
PPsvkn_2121 <- PP[SS=="Supm" & VV=="Batik" & KK=="High" & NN=="Low"]
PPsvkn_2122 <- PP[SS=="Supm" & VV=="Batik" & KK=="High" & NN=="High"]
PPsvkn_2211 <- PP[SS=="Supm" & VV=="Mariette" & KK=="Low" & NN=="Low"]
PPsvkn_2212 <- PP[SS=="Supm" & VV=="Mariette" & KK=="Low" & NN=="High"]
PPsvkn_2221 <- PP[SS=="Supm" & VV=="Mariette" & KK=="High" & NN=="Low"]
PPsvkn_2222 <- PP[SS=="Supm" & VV=="Mariette" & KK=="High" & NN=="High"]
svkn_ind[PPsvkn_1111] <- 16; svkn_ind[PPsvkn_1112] <- 15
svkn_ind[PPsvkn_1121] <- 14; svkn_ind[PPsvkn_1122] <- 13
svkn_ind[PPsvkn_1211] <- 12; svkn_ind[PPsvkn_1212] <- 11
svkn_ind[PPsvkn_1221] <- 10; svkn_ind[PPsvkn_1222] <- 9
svkn_ind[PPsvkn_2111] <- 8; svkn_ind[PPsvkn_2112] <- 7
svkn_ind[PPsvkn_2121] <- 6; svkn_ind[PPsvkn_2122] <- 5
svkn_ind[PPsvkn_2211] <- 4; svkn_ind[PPsvkn_2212] <- 3
svkn_ind[PPsvkn_2221] <- 2; svkn_ind[PPsvkn_2222] <- 1
X11()
par(font.axis=1,font.lab=1,font.main=1,pty="s")
plot(y=(svkn_ind[WW=="Fluct"]+0.1),x=disp_resp2[WW=="Fluct"],ylim=c(1,16),lab=c(5,16,7),
          yaxt="n",ylab="SVKN",cex.axis=1.2,cex.lab=1.4,xlim=c(0,0.5),cex=1.2,
          xlab="Squared Residuals")
axis(side=2,at=rev(c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)),
labels=c("1111","1112","1121","1122","1211","1212","1221","1222",
"2111","2112","2121","2122","2211","2212","2221","2222"),las=1)
points(y=(svkn_ind[WW=="Control"]-0.1),disp_resp2[WW=="Control"],pch=19,cex=1.2)


# -----------------------------------------------------------
# Figure 4
# -----------------------------------------------------------

# plot layout
X11()
nnf <- layout(matrix(c(1,2,3,
                      1,4,5,
                      0,6,6),3,3,byrow=TRUE),
                      c(1,5,5), c(5,5,1), TRUE)
layout.show(nnf)
# mean treatment effects
treat_VTR <- tapply(rate,list(VV,TT,RR),mean)
treat_111 <- treat_VTR["Batik","16C","300lux"]
treat_112 <- treat_VTR["Batik","16C","600lux"]
treat_121 <- treat_VTR["Batik","21C","300lux"]
treat_122 <- treat_VTR["Batik","21C","600lux"]
treat_211 <- treat_VTR["Mariette","16C","300lux"]
treat_212 <- treat_VTR["Mariette","16C","600lux"]
treat_221 <- treat_VTR["Mariette","21C","300lux"]
treat_222 <- treat_VTR["Mariette","21C","600lux"]
# y label
par(font.axis=1,font.lab=1,font.main=1,pty="m",mai=c(0,0,0,0),
       bty="n",xaxt="n",yaxt="n")
plot(x=1:2,y=1:2,pch=" ",ylim=c(1,2),xlim=c(1,2))
text(x=1.3,y=1.5,"Quality Score",srt=90,cex=2)
text(x=c(rep(1.85,4)),y=c(1.470,1.400,1.048,0.975)+0.008,
                 c(4,3,2,1),cex=1.8)
text(x=c(rep(1.85,4)),y=c(1.470,1.400,1.048,0.975)+0.550,
                 c(4,3,2,1),cex=1.8)
# plot 1
par(font.axis=1,font.lab=1,font.main=1,pty="s",mai=c(0,0,0,0),
       bty="o",xaxt="n",yaxt="n",tcl=0.3)
plot(x=seq(0,6,by=0.5),y=pred_funct2(treat_111,seq(-3,3,by=0.5))$prob,lwd=2,
         las=1,cex.axis=1.3,cex.lab=1.3,cex=1.4,xlab="Time (Years)",ylim=c(0,1),
             xlim=c(-0.5,6.5),type="b")
lines(x=seq(0,6,by=0.5),y=pred_funct2(treat_121,seq(-3,3,by=0.5))$prob,lwd=2,lty=1,
               type="b",pch=19,cex=1.4)
xtcks_fun(7)
yticks_fun(x=4,ytick=ord_y)
lines(x=c(-1,7),y=c(ord_y[2],ord_y[2]),lty=2,lwd=1,type="l")
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_111)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_121)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
text(x=6,y=1,"V1:R-",cex=1.8)
# plot 2
par(font.axis=1,font.lab=1,font.main=1,pty="s",mai=c(0,0,0,0),
       bty="o",xaxt="n",yaxt="n",tcl=0.3)
plot(x=seq(0,6,by=0.5),y=pred_funct2(treat_112,seq(-3,3,by=0.5))$prob,lwd=2,
         las=1,cex.axis=1.3,cex.lab=1.3,cex=1.3,xlab="Time (Years)",ylim=c(0,1),
             xlim=c(-0.5,6.5),type="b")
lines(x=seq(0,6,by=0.5),y=pred_funct2(treat_122,seq(-3,3,by=0.5))$prob,lwd=2,lty=1,
               type="b",pch=19,cex=1.4)
xtcks_fun(7)
yticks_fun(x=4,ytick=ord_y)
lines(x=c(-1,7),y=c(ord_y[2],ord_y[2]),lty=2,lwd=1,type="l")
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_112)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_122)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
text(x=6,y=1,"V1:R+",cex=1.8)
# plot 3
par(font.axis=1,font.lab=1,font.main=1,pty="s",mai=c(0,0,0,0),
       bty="o",xaxt="n",yaxt="n",tcl=0.3)
plot(x=seq(0,6,by=0.5),y=pred_funct2(treat_211,seq(-3,3,by=0.5))$prob,lwd=2,
         las=1,cex.axis=1.3,cex.lab=1.3,cex=1.3,xlab="Time (Years)",ylim=c(0,1),
             xlim=c(-0.5,6.5),type="b")
lines(x=seq(0,6,by=0.5),y=pred_funct2(treat_221,seq(-3,3,by=0.5))$prob,lwd=2,lty=1,
               type="b",pch=19,cex=1.4)
xtcks_fun(7)
yticks_fun(x=4,ytick=ord_y)
lines(x=c(-1,7),y=c(ord_y[2],ord_y[2]),lty=2,lwd=1,type="l")
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_211)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_221)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
text(x=6,y=1,"V2:R-",cex=1.8)
# plot 4
par(font.axis=1,font.lab=1,font.main=1,pty="s",mai=c(0,0,0,0),
       bty="o",xaxt="n",yaxt="n",tcl=0.3)
plot(x=seq(0,6,by=0.5),y=pred_funct2(treat_212,seq(-3,3,by=0.5))$prob,lwd=2,
         las=1,cex.axis=1.3,cex.lab=1.3,cex=1.3,xlab="Time (Years)",ylim=c(0,1),
             xlim=c(-0.5,6.5),type="b")
lines(x=seq(0,6,by=0.5),y=pred_funct2(treat_222,seq(-3,3,by=0.5))$prob,lwd=2,lty=1,
               type="b",pch=19,cex=1.4)
xtcks_fun(7)
yticks_fun(x=4,ytick=ord_y)
lines(x=c(-1,7),y=c(ord_y[2],ord_y[2]),lty=2,lwd=1,type="l")
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_212)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
lines(x=rep(ipred_funct2(prob=0.15,coeff=treat_222)$week,2),y=c(-0.05,ord_y[2]),lty=2,lwd=1)
text(x=6,y=1,"V2:R+",cex=1.8)
# x label
par(font.axis=1,font.lab=1,font.main=1,pty="m",mai=c(0,0,0,0),
       bty="n",xaxt="n",yaxt="n")
plot(x=1:2,y=1:2,pch=" ",ylim=c(1,2),xlim=c(1,2))
text(x=1.5,y=1.3,"Time in Home-life (weeks)",srt=0,cex=2)
text(x=c(1.020,1.090,1.160,1.230,1.305,1.375,1.445),
                y=rep(1.82,7),c(0,1,2,3,4,5,6),cex=1.8)
text(x=c(1.020,1.090,1.160,1.230,1.305,1.375,1.445)+0.539,
                y=rep(1.82,7),c(0,1,2,3,4,5,6),cex=1.8)


#------------------------------------------------------------
# End
#------------------------------------------------------------