varPredictNlmeGnls

Description

Predictions including variance for basic nlme and gnls models.

Usage

1
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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
#----  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=tmp <- 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