Topic: Growth

Topic Description:

Functions for analyzing tree growth, designed for the standard CTFS R Analytical Tables.

File: growth/growthfit.graph.r
View File Source Download File No help file available

Function: graph.growthmodel.spp

Function Description: graph.growthmodel.spp

Use output of growth.flexbin to graph observed growth and predictions. With conf>0, random draws from posterior parameters are used to make predictions and overlay with gray lines, and mean predicted value at each x is calculated. The option whichpred determines whether to graph the fitted values based on the median parameters of the Gibbs sampler, or the mean fitted value of all the steps in the Gibbs sampler.

Function Arguments:

ArgumentDefault Value
fit
jiggle.001
whichpred'pred'
xrangeNULL
yrangeNULL
xtitleNULL
ytitleNULL
includeaxsTRUE
'red'
regclr"green"
modelclr"blue"
modellwd1
graphdiv10
addFALSE
addptsTRUE
maintitleNULL
conf0

Function Source:


graph.growthmodel.spp=function(fit,jiggle=.001,whichpred='pred',xrange=NULL,yrange=NULL,xtitle=NULL,
        ytitle=NULL,includeaxs=TRUE,
withSD='red',regclr="green",modelclr="blue",modellwd=1,graphdiv=10,ad
        d=FALSE,addpts=TRUE,maintitle=NULL,conf=0) {
size=fit$model$x
obs=fit$model$y
if(whichpred=='meanpred') pred=fit$model$meanpred
else pred=fit$model$pred

if(is.null(yrange)) yrange=range(c(obs,pred))
if(is.null(xrange)) xrange=range(size)

if(is.null(xtitle)) xtitle='log size'
if(is.null(ytitle)) ytitle='growth'

xjiggle=rnorm(length(size),mean=0,sd=jiggle)
yjiggle=rnorm(length(size),mean=0,sd=jiggle)

if(addpts)
{

if(!add) plot(size+xjiggle,obs+yjiggle,xlim=xrange,ylim=yrange,pch=16,cex=.8,xlab=xtitle,ylab=yti
        tle,axes=includeaxs) else points(size+xjiggle,obs+yjiggle,pch=16,cex=.8)
}
else
{
if(add) lines(size,pred,col=modelclr,lwd=modellwd)

else plot(size,pred,type='l',xlim=xrange,ylim=yrange,xlab=xtitle,ylab=ytitle,col=modelclr,lwd=mod
        ellwd,axes=includeaxs) }

b=seq(min(size),max(size)+.0001,len=graphdiv-1)
sizecat=cut(size,breaks=b,right=FALSE)

if(conf>0)
{
allparam=fit$fullparam[fit$keep,]
noparam=dim(allparam)[1]
r=sample(1:noparam,conf)
# browser()
for(i in 1:conf)
{
onepred=linearmodel.bin(fit$model$x,param=allparam[r[i],])
lines(fit$model$x,onepred,col='gray10',lwd=0.5)
if(i==1) meanpred=onepred
else meanpred=meanpred+onepred
}

lines(fit$model$x,meanpred/conf,lty='dashed',lwd=1.5,col='purple')
}

meansize=tapply(size,sizecat,mean)
meangrowth=tapply(obs,sizecat,mean)

if(!is.null(regclr)) lines(meansize,meangrowth,col=regclr,lwd=modellwd)
if(addpts) lines(size,pred,col=modelclr,lwd=modellwd)

if(!is.null(withSD))
{
lines(size,pred+fit$model$predSD,col=withSD)
lines(size,pred-fit$model$predSD,col=withSD)
}

ypos=max(obs)-0.3*diff(range(obs))
xpos=min(size)+0.1*diff(range(size))
if(!is.null(maintitle)) text(xpos,ypos,maintitle,cex=1.7,pos=4)

return(data.frame(meansize,meangrowth))
}