[R] interactions and GAM

[Beaulaton Laurent]

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:laurent.beaulaton at bordeaux.cemagref.fr

http://www.laurent-beaulaton.fr/ 
http://www.bordeaux.cemagref.fr/rabx/