PLD<-inverts
Expq2<-lme(log(PLD)~I(log(temp)-2.8)+I((log(temp)-2.8)^2),
   random=~I(log(temp)-2.8)|species, data=PLD, method="ML", 
   control = list(niterEM=50))
#extract parameter estimates
sigmaf.inv<-solve(vcov(Expq2))
gamma.hat<-fixef(Expq2)
r<-length(gamma.hat)
species.list<-unique(PLD$species)
phi<- unlist(as.vector( attributes(summary(Expq2)$apVar)$Pars))
sigmar.inv<-solve(Expq2$apVar)
#obtain Cj
Cj<-NULL
for (i in 1:length(species.list)) {
   temp.dat<-PLD[PLD$species != species.list[i],]
   cur.model<-lme(log(PLD)~I(log(temp)-2.8)+I((log(temp)- 2.8)^2),
      random=~I(log(temp)-2.8)|species, data=temp.dat, method="ML", 
      control = list(niterEM=50))
   gamma.red<-fixef(cur.model)
   C0F<-1/(r+1)*(gamma.hat-gamma.red)%*%sigmaf.inv%*%(gamma.hat-gamma.red)
   phi.minusj<-unlist(as.vector(attributes(summary(cur.model)$apVar)$Pars))
   q<-length(phi.minusj)
   CjR<-1/q*(phi-phi.minusj)%*%sigmar.inv%*%(phi-phi.minusj)
   cur.Cj<-1/(r+q+1)*((r+1)*C0F+q*CjR)
   Cj<-c(Cj,cur.Cj)
}
#make covariance matrix of random effects
temp <- unlist(as.vector( attributes(summary(Expq2)$apVar)$Pars[1:3]))
D.hat <- matrix(0,2,2)
D.hat[1,1] <- exp(temp[1])^2
D.hat[2,2] <- exp(temp[2])^2
cor <-  2*(exp(temp[3])/(1+exp(temp[3])))-1
D.hat[1,2] <- cor*exp(temp[1])*exp(temp[2])
D.hat[2,1] <- D.hat[1,2]
sigusq.hat<-(as.numeric(exp(attributes(summary(Expq2)$apVar)$Pars[4])))^2
resids<-residuals(Expq2)
#calculate the generalized multivariate residuals
pval<-NULL
for (i in 1:length(species.list)) {
   cur.dat<-PLD[PLD$species == species.list[i],]
   cur.resids<-resids[names(resids)==species.list[i]]
   nj<-length(cur.dat$temp)
   xmat<-matrix(c(rep(1,nj),log(cur.dat$temp)-2.8),ncol=2)
   covmat<-xmat%*%D.hat%*%t(xmat)+diag(sigusq.hat,nrow=nj)
   inv.covmat<-solve(covmat)
   Sj2<-cur.resids%*%inv.covmat%*%cur.resids
   cur.pval<-1-pchisq(Sj2,nj)
   pval<-c(pval,cur.pval)
}
#plot results
legend.text<-substitute(expression(alpha == 0.05, 
   paste(alpha == x, ", Bonferroni")),
   list(x=round(.05/length(species.list),5)))
plot(pval,Cj,xlab="p-value of the standardized multivariate residuals",
   ylab=expression(paste("influence diagnostic  ",C[j],sep=' ')))
abline(v=.05/length(species.list),col=2,lty=3,lwd=2)
abline(v=.05,col=4,lty=2,lwd=2)
p.x<-round(pval[species.list=="Circeus spirillum"],4)
text(.115,.83,paste("p = ",p.x),cex=.7)
legend(.65,.85,legend=eval(legend.text),lty=c(2,3),lwd=c(2,2),col=c(4,2),
   bty='n',cex=c(.8,.8)) 
identify(pval,Cj,species.list,cex=.8,font=3)
mtext("Influence and Fit of Level-2 Units",side=3,
   line=.5,font=2,cex=.9)
rug(Cj,side=2,col='gray40',lwd=.05,ticksize=.02)
rug(pval,side=1,col='gray40',lwd=.05,ticksize=.02)
box()