Thanks a lot for all the mails!

With all the different examples I now have got it working. The problem
finally after I had worked out the newd2 as supplied by Gavin was that I
found the plot so strange looking. I have altered my newd2 a bit, and after
your e-mails Simon I figured out how to finally plot my different graphs.
The final script I used to get things running is as follows, maybe a bit
different then the examples but I think it did give me the predictions I
wanted. This is the version I used for the plots for all days.


m3<-gam(y_no~s(day,by=mapID)+mapID,family=binomial, data=mergeday)

newd2 <- data.frame(day = rep(seq(1, 366, length = 366), 8),
                   mapID = rep(sort(unique(mergeday$mapID)),
                               each = 366))

pred<-data.frame(predict(m3,newd2,type="response"))

names(pred)
[1] "predict.m3..newd2..type....response.."

### I need to rename the column name in something more sensible or easy, but
this worked out as well.
plot(pred$predict.m3..newd2..type....response..[1:366],type="l",ylab="Relative
Habitat Use",xlab ="days")

If this method might not give me a correct prediction, please tell me, but I
think the results looked kind of similar to my GAM plots.

A mighty thank you for all the helpful replies!

Robbert


2009/1/21 Simon Wood <s.wood@bath.ac.uk>

> 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:
> > Thanks for the swift reply,
> >
> > I might have been a bit sloppy with describing my datasets and problem. I
> > showed the first model as an example of the type of GAM that I had been
> > able to use the predict function on. What I am looking for is how to
> > predict my m3:
> > model3<-gam(y_no~s(day,by=mapID),family=binomial, data=mergeday)
> >
> > When I plot this I get 8 different graphs. Each showing me a different
> > habitat type with on the x-axis days and on the y-axis s(day,2,81):mapID.
> > With predict I was hoping to get the scale of the y-axis right for a
> > selection of days (for example 244,304).
> >
> > I have tried to reform the script you gave me to match my dataset in m3,
> > but it all did not seem to work.
> >
> > newd2 <- data.frame(day = rep(seq(244, 304, length = 100), 8),
> >                    mapID = rep(levels(mergeday$mapID), each = 100))
> >
> > newd2 <- data.frame(day = rep(seq(244, 304, length = 100), 8),
> >                    mapID = rep(sort(unique(mergeday$mapID)),
> >                                each = 100))
> >
> > I am guessing it must have something to do with the  by in
> s(day,by=mapID).
> > I haven't come across any examples that used a GAM with by and then used
> > the predict function.
> >
> > (A sample of the dataset:
> >           mapID                  day   y_no
> > Urban Areas and Water    25      1
> > Urban Areas and Water    26      1
> > Early Succesional Forest  27      0
> > Agriculture                       28      0
> > Early Succesional Forest  29      0
> > Mature Coniferous Forest  30      0)
> >
> >
> > I am sorry that I have to bother you even more with this, and I hope that
> > my additional explanation about my problem might help solve it.
> >
> > Sincerely yours,
> >
> > Robbert Langenberg
> >
> > 2009/1/16 Gavin Simpson <gavin.simpson@ucl.ac.uk>
> >
> > > On Fri, 2009-01-16 at 12:36 +0100, Robbert Langenberg wrote:
> > > > Dear,
> > > >
> > > > I am trying to get a prediction of my GAM on a response type. So that
> I
> > > > eventually get plots with the correct values on my ylab.
> > > > I have been able to get some of my GAM's working with the example
> shown
> > > > below:
> > > > *
> > > > model1<-gam(nsdall ~ s(jdaylitr2), data=datansd)
> > > > newd1 <- data.frame(jdaylitr2=(244:304))
> > > > pred1 <- predict.gam(model1,newd1,type="response")*
> > >
> > > 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
> > >
> > > > The problem I am encountering now is that I cannot seem to get it
> done
> > >
> > > for
> > >
> > > > the following type of model:
> > > >
> > > > *model3<-gam(y_no~s(day,by=mapID),family=binomial, data=mergeday)*
> > > >
> > > > My mapID consists of 8 levels of which I get individual plots with *
> > > > plot(model3)*. When I do predict with a newdata in it just like my
> > > > first model I need all columns to have the same amount of rows or
> else
> > > > R will
> > >
> > > not
> > >
> > > > except it ofcourse, the col.names need to at least include day and
> > > > mapID. This way I can not get a prediction working for this GAM, I am
> > > > confused because of this part in the model: *s(day,by=mapID).
> > > >
> > > > *I have been reading through the GAM, an introduction with R book
> from
> > >
> > > Wood,
> > >
> > > > S. but could not find anything about predictions with BY in the
> model.
> > > >
> > > > I hope someone can help me out with this,
> > > >
> > > > Sincerely yours,
> > > >
> > > > Robbert Langenberg
> > > >
> > > >       [[alternative HTML version deleted]]
> > > >
> > > > ______________________________________________
> > > > R-help@r-project.org 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.
> > >
> > > --
> > > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
> > >  Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
> > >  ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
> > >  Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
> > >  Gower Street, London          [w]
> > > http://www.ucl.ac.uk/~ucfagls/ <http://www.ucl.ac.uk/%7Eucfagls/><
> http://www.ucl.ac.uk/%7Eucfagls/> UK. WC1E
> > > 6BT.                 [w] http://www.freshwaters.org.uk
> > > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
>
> --
> > Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> > +44 1225 386603  www.maths.bath.ac.uk/~sw283<http://www.maths.bath.ac.uk/%7Esw283>
>



-- 
www.lowlandpaddies.nl

	[[alternative HTML version deleted]]


