# R-script that produces figure 4 of paper 
# Mehtätalo, L., Peltola, H., Kilpeläinen, A. and Ikonen, V.-P. 2014. The response of 
# basal area growth of Scots pine to thinning: A longitudinal analysis of tree-specific 
# series using a nonlinear mixed-effects model. Forest Science 60 (4): pp. 636-644. 
# DOI: http://dx.doi.org/10.5849/forsci.13-059. 
# 
# Lauri Mehtätalo 3.3.2017

load("c:/laurim/patti/fig4.Rdata")
library(nlme)
ap<-augPred(ThinNlme91)

trees<-c(108,408,808,2610,11910,23710,11305,17605,23905,
		106,406,306,706,2104,6204,9404,13404,
		103,303,603,803,807,2807,4807,8907)[c(7,9,2,4,8,5,6,1,3,16,11,15,12,14,9,13,17,18,22,19,24,23,21,20,25)]

patti4$pred0<-predict(ThinNlme91,level=0)
patti4$pred<-predict(ThinNlme91)

# pdf("d:/laurim/patti/Fig4.pdf")
treatments<-c("Control","Light","Moderate","Heavy")
windows()
layout(matrix(1:36, 6, 6, byrow = TRUE),widths=c(1,2.1,2.1,2.1,2.1,2.1),heights=c(1.5,1.5,1.5,1.5,1.5,0.6))
for (i in 1:length(trees)) {
	if (i==1|i==6|i==11|i==16|i==21) {
		par(mai=c(0,0.0,0,0),col.axis="white",mgp=c(1.6,0.8,0),bty="n",xaxt="n",yaxt="n")
		plot(1,type="n",ylim=range(patti4$ThEf),xlim=c(-1,3)) 
		text(1.8,c(-200,0,200,600),c(-200,0,200,600),srt=90)
		if (i==11) text(0.6,200,"ThEf, mm2/yr",srt=90,cex=1.2)
	} 
	this<-patti4[patti4$Tree==trees[i],]
	par(mai=c(0,0.0,0.05,0),col.axis="white",mgp=c(1.6,0.8,0),bty="l",xaxt="s",yaxt="s",ann=FALSE,cex.axis=0.01) 
	plot(this$ythin,this$ThEf,cex=0.1,ylim=range(patti4$ThEf),type="l",col=gray(0.5))
	text(10,600,paste(treatments[this$SDClass[1]]))
	text(10,500,paste("DBH=",round(this$Diam1986[1]/10,1)))
	lines(this$ythin,this$pred0,lty="dashed")
	lines(this$ythin,this$pred)
}
par(mai=c(0,0.0,0,0),col.axis="white",mgp=c(1.6,0.8,0),bty="n",xaxt="n",yaxt="n")
plot(1,type="n",xlim=c(0,18),ylim=c(-1,3)) 
for (i in 1:5) {
	plot(1,type="n",xlim=c(0,18),ylim=c(-1,3)) 
	text(c(0,5,10,15),2,c(0,5,10,15))
	if (i==3) text(9,0.5,"Time since thinning",cex=1.2)
}
#dev.off()