Don't know if this helps, but... gam in package mgcv will let you set up smooths that interact with factors using the `by' variable mechanism. See ?gam.models, particularly the last example. Prediction is supported.
Simon -- > Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK > +44 1225 386603 www.maths.bath.ac.uk/~sw283 On Tuesday 27 February 2007 20:10, Beaulaton Laurent wrote: > Dear R-users, > > I have 1 remark and 1 question on the inclusion of interactions in the gam > function from the gam package. > > I need to fit quantitative predictors in interactions with factors. You can > see an example of what I need in fig 9.13 p265 from Hastie and Tibshirani > book (1990). It's clearly stated that in ?gam "Interactions with > nonparametric smooth terms are not fully supported". I have found a trick > in a former http://www.math.yorku.ca/Who/Faculty/Monette/S-news/2284.html, > using NAs and na.gam.replace argument, but some points are still unclear > for me. > > First the prediction of new data (using predict function) is not so easy > (see script below), and need a close reading from section 7.3.2 of the > Chambers and Hastie (1992). > > Second I need to have the same intercept for all levels of factor and this > not achievable with this trick. My question is : why not replacing NA by 0 > (or any other particular value) ? > > Here is a quite long (sorry for that) script with a generated dataset to > better undestand my question. in this script the model to fit is (in a > GLM-like writing) : y~s(x2):x1 the generated dataset follows this model and > y(x2=0)=10 whatever x1. > > ######################## > #start of script > ######################## > > #data construction (with deliberately very small noise) > data1=data.frame(x1=rep(NA,27),x2=NA,y=NA) > > data1$x1=factor(c(rep(1,11),rep(2,11),rep(3,5))) > data1$x2=c(rep(0:10,2),0:4) > > data1[data1$x1==1,"y"]=data1[data1$x1==1,"x2"]^4*5+rnorm(11)+10000 > data1[data1$x1==2,"y"]=data1[data1$x1==2,"x2"]^4*(-3)+rnorm(11)+10000 > data1[data1$x1==3,"y"]=10000*data1[data1$x1==3,"x2"]+rnorm(5)+10000 > > library(lattice) > xyplot(data1$y~data1$x2,groups=data1$x1) > > #creation of dummy variables for interactions > data1$x2_1=ifelse(data1$x1=="1",data1$x2,NA) > data1$x2_2=ifelse(data1$x1=="2",data1$x2,NA) > data1$x2_3=ifelse(data1$x1=="3",data1$x2,NA) > > #model fitting > library(gam) > model=gam(y~s(x2_1)+s(x2_2)+s(x2_3)+x1,data=data1,na=na.gam.replace) > > #prediction fit well data : > summary(model) > plot(data1$x2,data1$y) > points(data1$x2,model$fitted.value,col="red",pch="+") > > #trying to see prediction: > predict(model) #does work > predict(model,newdata=data1) #produce NA > > #trying to replace NA in data1 by mean, to mimic na.gam.replace: > Ndata=data1 > Ndata$x2_1=ifelse(data1$x1=="1",data1$x2,mean(data1$x2_1,na.rm=TRUE)) > Ndata$x2_2=ifelse(data1$x1=="2",data1$x2,mean(data1$x2_2,na.rm=TRUE)) > Ndata$x2_3=ifelse(data1$x1=="3",data1$x2,mean(data1$x2_3,na.rm=TRUE)) > > predict(model,Ndata)-predict(model) #as you can see there is a systematic > biais > > #correct way to predict (=returned 0 for terms with NA value): > p=predict(model,data1,type="term") > rowSums(cbind(p,attr(p,"constant")),na.rm=TRUE)-predict(model) > > #alternative solution, 0 instead of NA > data1$v1=ifelse(data1$x1=="1",data1$x2,0) > data1$v2=ifelse(data1$x1=="2",data1$x2,0) > data1$v3=ifelse(data1$x1=="3",data1$x2,0) > > model1=gam(y~s(v1)+s(v2)+s(v3),data=data1) > summary(model1) > points(data1$x2,predict(model1,data1),col="green",pch="X") > #no particular problem with predict function > > #what's happened in x2=0 ? > predict(model)[data1$x2==0] > predict(model1)[data1$x2==0] > > ######################## > #end of script > ######################## > > thanks in advance > best regards > Laurent Beaulaton > > --------------------------------------------- > Laurent Beaulaton > ############################### > # NEW !!!! # > # http://www.laurent-beaulaton.fr/ # > # Tel + 33 (0)5 57 89 27 17 # > ############################### > --------------------------------------------- > Cemagref (French Institute of Agricultural and Environmental Engineering > Research ) Unité "Ecosystèmes estuariens et poissons migrateurs > amphihalins" > (anciennement Unité "Ressources aquatiques continentales") > 50 avenue de Verdun > F 33612 Cestas Cedex > > Tel + 33 (0)5 57 89 27 17 > Fax + 33 (0)5 57 89 08 01 > mailto:[EMAIL PROTECTED] > > http://www.laurent-beaulaton.fr/ > http://www.bordeaux.cemagref.fr/rabx/ > > ______________________________________________ > [email protected] 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. ______________________________________________ [email protected] 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.
