In addition to this, is it in anyway possible to get the variance in the prediction graph as a shade similar to the GAM plots?
2009/1/21 Robbert Langenberg <mcre...@gmail.com> > 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 > -- 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.
