[R] gam parameter predictions --Sorry for double posting

Simon Wood s.wood at bath.ac.uk
Tue Mar 27 13:29:46 CEST 2007


Looks ok to me, provided that you want averages (on log scale) taken over the 
the observed covariate values. If you want variances for the means you could 
always do something like the following...

Xa <- model.matrix(~as.factor(year)-1)
Xa <- t(Xa)/colSums(Xa)  ## Xa%*%fitted(g1) gives required averages

Xp <- predict(g1,type="lpmatrix") ## Xp%*%coef(g1) gives fitted values
Xm <- Xa%*%Xp  ## now Xm%*%coef(g1) gives required averages

yearly.means <- Xm%*%coef(g1)
ym.cov <- Xm%*%vcov(g1,freq=FALSE)%*%t(Xm) ## cov matrix of yearly means

(help file for predict.gam has a bit more information on this sort of thing, 
or section 5.2.6 of my book)

best,
Simon

-- 
> Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> +44 1225 386603  www.maths.bath.ac.uk/~sw283 

>
> The covariates are year, month, vessel and statistical rectangle.
>
>
> The model looks like this:
>
> g1 <- gam(log(cpue) ~  s(rekt1) + s(year) + s(mon) + s(reg1), data =
> dataTest)
>
>
> Once the model is fitted to the data I want to get the mean model
> estimates by year.
>
> I do the following:
>
> obsPred <- data.frame(year = dataTest$year, pred = predict(g1, type =
> "response"))
>
> gamFit <- tapply(obsPred$pred, list(year = obsPred$ar), mean)
>
>
>
> Is this correct?
>
>
>
> Thanks in advance
>
> > version
>
>                _
> platform       i386-pc-mingw32
> arch           i386
> os             mingw32
> system         i386, mingw32
> status
> major          2
> minor          4.1
> year           2006
> month          12
> day            18
> svn rev        40228
> language       R
> version.string R version 2.4.1 (2006-12-18)
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented, minimal,
> self-contained, reproducible code.



More information about the R-help mailing list