resids.table<-function(DATASET,Model.out)
{
#get list of species names
species.list<-unique(DATASET$species)
#construct the matrix Omega
tau00<-as.numeric((VarCorr(Model.out))[1,1])
tau11<-as.numeric((VarCorr(Model.out))[2,1])
rho01<-as.numeric((VarCorr(Model.out))[2,3])
tau01<-sqrt(tau00*tau11)*rho01
omegaH<-matrix(c(tau00,tau01,tau01,tau11),ncol=2)
#obtain the level-2 residuals from the model
reslev2.u0<-ranef(Model.out)[,1] 
reslev2.u1<-ranef(Model.out)[,2]

#limma is needed for the blockDiag function
library(limma)

#locate where each level-2 unit starts and stops
cumsum(table(DATASET$species))->stops
starts<-c(1,stops+1)
starts<-starts[-length(starts)]

#initialize the matrices V, R, and S
Vmat<-NULL
Rmat<-NULL
sh<-NULL

#build the block diagonal matrices species by species
for (i in 1:length(stops)) {
  #next line creates Zh--currently an intercept and lntemp-2.8
   Zmat<-cbind(rep(1,stops[i]-starts[i]+1),
       DATASET[starts[i]:stops[i],"lntemp"]-2.8) 
   Rhat<-Zmat%*%omegaH
   Vpart<-Zmat%*%omegaH%*%t(Zmat)  
   Vmat<-if(i==1) Vpart else blockDiag(Vmat,Vpart)
   Rmat<-if(i==1) Rhat else blockDiag(Rmat,Rhat)
   sh<-if(i==1) omegaH else blockDiag(sh,omegaH)
  }   
#add sigma2 to the diagonal entries
sigma2<-as.numeric((VarCorr(Model.out))[3,1])
V.true<-Vmat+diag(rep(sigma2,dim(Vmat)[1]))

#need X matrix of model variables. Currently these are
#intercept, lntemp-2.8, and (lntemp-2.8)^2
xmat<-cbind(rep(1,dim(DATASET)[1]),
   DATASET[,"lntemp"]-2.8,(DATASET[,"lntemp"]-2.8)^2)

#sand.wich is the middle part of variance expression—the
#part in the parentheses plus premultiplication by V^-1
sand.wich<-(solve(V.true))%*%(V.true-xmat%*%
   (solve(t(xmat)%*%(solve(V.true))%*%xmat))%*%t(xmat))
res2.sd<-sqrt(diag(sh-t(Rmat)%*%sand.wich%*%(solve(V.true))%*%Rmat))

#pull out the different random effects
sd.u0<-res2.sd[names(res2.sd)=='1']
sd.u1<-res2.sd[names(res2.sd)=='2']

#construct data frame of results
out.resids<-data.frame(as.numeric(species.list),reslev2.u0,
   reslev2.u1,sd.u0,sd.u1)
colnames(out.resids)<-c("Species.num","u0","u1","sd.u0","sd.u1")
rownames(out.resids)<-species.list
out.resids
}
#### end function ####

out.resids<-resids.table(inverts,model2)



graph.u1<-function(out.resids,dataset) {
u1.upper<-out.resids$u1+qnorm(.975)*out.resids$sd.u1
u1.lower<-out.resids$u1+qnorm(.025)*out.resids$sd.u1
species.list<-unique(dataset$species)
plot(out.resids[order(out.resids$u1),"u1"],pch=17,cex=.8,
   ylim=c(min(u1.lower),max(u1.upper)),
   xlab='Rank order of Residuals',
   ylab=expression(paste('Level-2 Residual  ',u["1i"])))
sort.lower<-u1.lower[order(out.resids$u1)]
sort.upper<-u1.upper[order(out.resids$u1)]
for(i in 1:(dim(out.resids)[1]))
arrows(i,sort.lower[i],i,sort.upper[i],
  angle=90,code=3,length=.05,col='gray60')
abline(h=0,col=4,lty=2,lwd=2)
points(1:length(species.list),out.resids[order(out.resids$u1),"u1"],
   pch=17,cex=.8)
identify(1:length(species.list),out.resids[order(out.resids$u1),"u1"],
   rownames(out.resids[order(out.resids$u1),]),cex=.7,font=4,col=2)
mtext(expression(paste(bold('Caterpillar Plot for  '),bold(u["1i"]))),
   side=3,line=.5,font=2)
}
graph.u1(out.resids,inverts)



graph.u0<-function(out.resids,dataset) {
u0.upper<-out.resids$u0+qnorm(.975)*out.resids$sd.u0
u0.lower<-out.resids$u0+qnorm(.025)*out.resids$sd.u0
species.list<-unique(dataset$species)
plot(out.resids[order(out.resids$u0),"u0"],pch=17,cex=.8,
   ylim=c(min(u0.lower),max(u0.upper)),
   xlab='Rank order of Residuals',
   ylab=expression(paste('Level-2 Residual  ',u["0i"])))
sort.lower<-u0.lower[order(out.resids$u0)]
sort.upper<-u0.upper[order(out.resids$u0)]
for(i in 1:(dim(out.resids)[1]))
arrows(i,sort.lower[i],i,sort.upper[i],
   angle=90,code=3,length=.05,col='gray60')
abline(h=0,col=4,lty=2,lwd=2)
points(1:length(species.list),out.resids[order(out.resids$u0),"u0"],
   pch=17,cex=.8)
identify(1:length(species.list),out.resids[order(out.resids$u0),"u0"],
   rownames(out.resids[order(out.resids$u0),]),cex=.7,font=3,col=2)
mtext(expression(paste(bold('Caterpillar Plot for  '),bold(u["0i"]))),
   side=3,line=.5,font=2)
}
graph.u0(out.resids,inverts)