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.w...@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.simp...@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]]
______________________________________________
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.
