Skip to content

Predictions with GAM

9 messages · Gavin Simpson, Ken Knoblauch, Simon Wood +1 more

Original
Gavin Simpson
#
On Fri, 2009-01-16 at 12:36 +0100, Robbert Langenberg wrote:
Hi Robert,

You want predictions for the covariate over range 244:304 for each of
your 8 mapID's, yes?

This is not tested, but why not something like:

newd2 <- data.frame(jdaylitr2 = rep(seq(244, 304, length = 100), 8),
                    mapID = rep(levels(datansd$mapID), each = 100))

Then use newd2 in your call to predict.

I am assuming that datansd$mapID is a factor in the above. If it is just
some other indicator variable, then perhaps something like:

newd2 <- data.frame(jdaylitr2 = rep(seq(244, 304, length = 100), 8),
                    mapID = rep(sort(unique(datansd$mapID)), 
                                each = 100))

Does that work for you?

HTH

G
Ken Knoblauch
# 
Hi,

Robbert Langenberg <mcrelay <at> gmail.com> writes:
I'm not sure that this will work for you, but I had a similar
situation and was able to get predict to work (after helpful
advice from Simon Wood) with a by variable by generating 
a model matrix for a model with the interaction of the covariate 
and the by term, something like

model.matrix(~ day:mapID - 1, data = mergeday)

in your case.
I added the appropriate columns into my data frame and 
also to the newdata for predict.  You can see an example 
in the appendix of

http://www.journalofvision.org/8/16/10/

HTH,

Ken
Ken Knoblauch
# 
Ken Knoblauch <ken.knoblauch <at> inserm.fr> writes:
Let me correct this before someone else does for me,
(It's getting late on Friday afternoon!)
The interaction in the model matrix would not be
with day.  In my case, I had a covariate for the by
variable and a factor with two levels, so I used
the interaction to separate them with respect to
the smooths.  I think that you might be able to get
the columns from the model matrix columns from
the mapID term by itself.  Sorry for the misstep.

Ken
2 days later
Simon Wood
# 
Robbert,
On Friday 16 January 2009 14:30, Robbert Langenberg wrote:
--- Can you explain how it didn't work? What exactly was wrong with what 
predict.gam gave you? [when I try predict.gam for a gam with factor by 
variables, I seem to get what I'd expect, so I need a bit more detail to 
figure out what the issue is.] 

-- incidentally, is there any reason for omitting the main effect of `mapID' 
from the model? I'd usually expect to see 
  y_no~s(day,by=mapID)+mapID
as the model formula in this sort of situation. (see ?gam.models for details)

best,
Simon

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      Simon Wood
    

    

      
      #
    

  


  

      
Here's a simulated example (although really for this model structure, one 
might as well  fit seperate models for each factor level). 

## Simulate some data with factor dependent smooths
n <- 400
x <- runif(n, 0, 1)
f1 <- 2 * sin(pi * x)
f2 <- exp(2 * x) - 3.75887
f3 <- 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 *
            (1 - x)^10
fac <- as.factor(c(rep(1, 100), rep(2, 100), rep(3, 200)))
fac.1 <- as.numeric(fac == 1)
fac.2 <- as.numeric(fac == 2)
fac.3 <- as.numeric(fac == 3)
f <- f1 * fac.1 + f2 * fac.2 + f3 * fac.3
y <- rpois(f,exp(f/4))

## fit gam, with a smooth of `x' for each level of `fac'
b <- gam(y~s(x,by=fac)+fac,family=poisson)
par(mfrow=c(2,2))
plot(b)

## produce plots on response scale, first the prediction...
np <- 200
newd <- data.frame(x=rep(seq(0,1,length=np),3),
                   fac=factor(c(rep(1,np),rep(2,np),rep(3,np))))
pv <- predict(b,newd,type="response")

## .. now the plotting
par(mfrow=c(2,2))
ind <- 1:np
plot(newd$x[ind],pv[ind],type="l",xlab="x",ylab="f(x,fac=1)")
ind <- ind+np
plot(newd$x[ind],pv[ind],type="l",xlab="x",ylab="f(x,fac=2)")
ind <- ind+np
plot(newd$x[ind],pv[ind],type="l",xlab="x",ylab="f(x,fac=2)")

        
On Friday 16 January 2009 14:30, Robbert Langenberg wrote: