varPredictNlmeGnls: varPredictNlmeGnls

In bgctw/twNlme: Standard errors for basic nlme and gnls model predictions

Description

Predictions including variance for basic nlme and gnls models.

Usage

varPredictNlmeGnls(object, newdata, ...)

Arguments

object	the model fit object used for predictions, treated by attachVarPrep
newdata	dataframe of new predictors and covariates
...	further arguments to predict.lme or predict.gls

Details

Variance calculation is based on Taylor series expansion as described in appendix A1 by Wutzler08.

Fitted object needs to be prepared by function attachVarPrep. If not done before, this function is called automatically within varPredictNlmeGnls. However, for finetuning or avoiding overhead in repeated calls, it is recommended to explicitely call attachVarPrep before calling varPredictNlmeGnls.

Value

numeric matrix with columns

fit	predictions
varFix	variance component due to uncertainty in fixed effects
varRan	variance component due to uncertainty in random effects
varResid	variance component due to residual error
sdPop	standard deviation of prediction of a new population
sdInd	standard deviation of prediction of a new individual

Author(s)

Thomas Wutzler

References

Wutzler, T.; Wirth, C. & Schumacher, J. (2008) Generic biomass functions for Common beech (Fagus sylvatica L.) in Central Europe - predictions and components of uncertainty. Canadian Journal of Forest Research, 38, 1661-1675

See Also

varSumPredictNlmeGnls, twNlme-package

Examples

#----  fit a nlme and gnls model to data of stem weights
#data(Wutzler08BeechStem)

lmStart <- lm(log(stem) ~ log(dbh) + log(height), Wutzler08BeechStem )
nlmeFit <- nlme( stem~b0*dbh^b1*height^b2, data=Wutzler08BeechStem
			,fixed=list(b0 ~ si + age + alt, b1+b2 ~ 1)
			,random=  b0 ~ 1 | author
			,start=c( b0=c(as.numeric(exp(coef(lmStart)[1])),0,0,0)
			          , b1=as.numeric(coef(lmStart)[2])
			          , b2=as.numeric(coef(lmStart)[3]) )
			,weights=varPower(form=~fitted(.))
			,method='REML'		# for unbiased error estimates
		)
summary(nlmeFit)
x3 <- update(nlmeFit, fixed=list(b0 ~ si * log(age), b1+b2 ~ 1))

gnlsFit <- gnls( stem~b0*dbh^b1*height^b2, data=Wutzler08BeechStem
	,params = list(b0 ~ si + age + alt, b1~1, b2 ~ 1)
	,start=c( b0=c(as.numeric(exp(coef(lmStart)[1])),0,0,0)
	          , b1=as.numeric(coef(lmStart)[2])
	          , b2=as.numeric(coef(lmStart)[3]) )
	,weights=varPower(form=~fitted(.))
)
summary(gnlsFit)
fixef(gnlsFit)		# note the usage of fixef.gnls.
ranef(gnlsFit)		# note the usage of ranef.gnls.

#---- some artificial data for new prediction
nData <- data.frame(
  dbh=seq(2,80,length.out=40)
  , alt=median(Wutzler08BeechStem$alt)
  , si=median(Wutzler08BeechStem$si) )
lmHeight <- lm( height ~ dbh, Wutzler08BeechStem)
nData$height <- predict(lmHeight, nData)
lmAge <- lm( age ~ dbh, Wutzler08BeechStem)
nData$age <- predict(lmAge, nData)

#---- do the prediction including variance calculation
# automatic derivation with accounting for residual variance model
nlmeFit <- attachVarPrep(nlmeFit, fVarResidual=varResidPower)
resNlme <- varPredictNlmeGnls(nlmeFit,nData)

# plotting prediction and standard errors
plot( resNlme[,"fit"] ~ dbh, nData, type="l", xlim=c(40,80), lty="dashed")
lines( resNlme[,"fit"]+resNlme[,"sdPop"] ~ dbh, nData, col="maroon", lty="dashed" )
lines( resNlme[,"fit"]-resNlme[,"sdPop"] ~ dbh, nData, col="maroon", lty="dashed" )
lines( resNlme[,"fit"]+resNlme[,"sdInd"] ~ dbh, nData, col="orange", lty="dashed"  )
lines( resNlme[,"fit"]-resNlme[,"sdInd"] ~ dbh, nData, col="orange", lty="dashed" )

#---- handling special model of residual weights
# here we fit different power coefficients for authors
# for the prediction we take the mean, but because it appears in a nonlinear term
# we need a correction term
.tmp.f <- function(){	# takes long, so do not execute each test time
	nlmeFitAuthor <- nlme( stem~b0*dbh^b1*height^b2, data=Wutzler08BeechStem
		,fixed=list(b0 ~ si + age + alt, b1+b2 ~ 1)
		,random=  b0 ~ 1 | author
		,start=c( b0=c(as.numeric(exp(coef(lmStart)[1])),0,0,0)
		          , b1=as.numeric(coef(lmStart)[2])
		          , b2=as.numeric(coef(lmStart)[3]) )
		,weights=varPower(form=~fitted(.)|author)
	)
	#pred <- predict(modExampleStem, newdata, level=0)
	varResidPowerAuthor <- function(object,newdata,pred	){
		sigma <- object$sigma
		deltaAuthor <- coef(object$modelStruct$varStruct, allCoef = TRUE)
		delta <-  mean(deltaAuthor)
		varDelta <- var(deltaAuthor)
		sigma^2 * abs(pred)^(2*delta) * (1+2*log(abs(pred))^2*varDelta)
	}
	#mtrace(varResidPowerAuthor)
	nlmeFitResidAuthor <- attachVarPrep(
	  nlmeFitAuthor, fVarResidual=varResidPowerAuthor)
	resNlmeAuthor <- varPredictNlmeGnls(nlmeFitResidAuthor,nData)
	AIC(nlmeFitResidAuthor)
}

nlmeFit	#return for creating modExampleStem.RData

bgctw/twNlme documentation built on May 30, 2019, 4:01 p.m.