hi! I hope all of you can help me this question for example GAMs: ozonea<-gam(newozone~ pressure+maxtemp+s(avetemp,bs="cr")+s(ratio,bs="cr"),family=gaussian (link=log),groupA,methods=REML) formula(ozonea) newozone ~ pressure + maxtemp + s(avetemp, bs = "cr") + s(ratio,bs = "cr") #formula of gams coef(ozonea) # extract the coefficient of GAMs (Intercept) pressure maxtemp s(avetemp).1 s(avetemp).2 s(avetemp).3 s(avetemp).4 s(avetemp).5 s(avetemp).6 s(avetemp).7 s(avetemp).8 s(avetemp).9 s(ratio).1 2.848204561 -0.024662520 -0.095528572 -0.052613005 -0.135309283 0.066243642 0.030937179 0.176267227 0.243405931 0.357027920 0.600628869 0.758581430 0.004653696 s(ratio).2 s(ratio).3 s(ratio).4 s(ratio).5 s(ratio).6 s(ratio).7 s(ratio).8 s(ratio).9 0.156547789 -0.132527837 -0.179947367 -0.208918785 -0.468014803 -0.195498523 -0.813667830 -1.505844690 then, I want to calculate newozone - s(ratio,bs="cr") we define this term (newozone - s(ratio,bs="cr")) is X can i use the coefficient to calculate? -- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3900848.html Sent from the R help mailing list archive at Nabble.com.
Question about GAMs
8 messages · Daniel Malter, Simon Wood, Gavin Simpson +2 more
I am not an expert on this, but there is a way to check this. You can predict from a gam using predict(ozonea, newdata=...). In the "newdata" argument you can specify the X-values of interest to you. Thus, you can compare if your predictions are the same when predicted directly from the gam or when predicted by hand. HTH, Daniel
pigpigmeow wrote:
hi! I hope all of you can help me this question for example GAMs: ozonea<-gam(newozone~ pressure+maxtemp+s(avetemp,bs="cr")+s(ratio,bs="cr"),family=gaussian (link=log),groupA,methods=REML) formula(ozonea) newozone ~ pressure + maxtemp + s(avetemp, bs = "cr") + s(ratio,bs = "cr") #formula of gams coef(ozonea) # extract the coefficient of GAMs (Intercept) pressure maxtemp s(avetemp).1 s(avetemp).2 s(avetemp).3 s(avetemp).4 s(avetemp).5 s(avetemp).6 s(avetemp).7 s(avetemp).8 s(avetemp).9 s(ratio).1 2.848204561 -0.024662520 -0.095528572 -0.052613005 -0.135309283 0.066243642 0.030937179 0.176267227 0.243405931 0.357027920 0.600628869 0.758581430 0.004653696 s(ratio).2 s(ratio).3 s(ratio).4 s(ratio).5 s(ratio).6 s(ratio).7 s(ratio).8 s(ratio).9 0.156547789 -0.132527837 -0.179947367 -0.208918785 -0.468014803 -0.195498523 -0.813667830 -1.505844690 then, I want to calculate newozone - s(ratio,bs="cr") we define this term (newozone - s(ratio,bs="cr")) is X can i use the coefficient to calculate?
-- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3901383.html Sent from the R help mailing list archive at Nabble.com.
I'd be inclined to use predict(ozonea,type="terms") to extract the estimates of s(ratio,bs = "cr") that you need. But do you really want newozone - s(ratio,bs="cr") when you've used a log link? best, Simon
On 10/13/2011 09:05 AM, pigpigmeow wrote:
hi! I hope all of you can help me this question for example GAMs: ozonea<-gam(newozone~ pressure+maxtemp+s(avetemp,bs="cr")+s(ratio,bs="cr"),family=gaussian (link=log),groupA,methods=REML) formula(ozonea) newozone ~ pressure + maxtemp + s(avetemp, bs = "cr") + s(ratio,bs = "cr") #formula of gams coef(ozonea) # extract the coefficient of GAMs (Intercept) pressure maxtemp s(avetemp).1 s(avetemp).2 s(avetemp).3 s(avetemp).4 s(avetemp).5 s(avetemp).6 s(avetemp).7 s(avetemp).8 s(avetemp).9 s(ratio).1 2.848204561 -0.024662520 -0.095528572 -0.052613005 -0.135309283 0.066243642 0.030937179 0.176267227 0.243405931 0.357027920 0.600628869 0.758581430 0.004653696 s(ratio).2 s(ratio).3 s(ratio).4 s(ratio).5 s(ratio).6 s(ratio).7 s(ratio).8 s(ratio).9 0.156547789 -0.132527837 -0.179947367 -0.208918785 -0.468014803 -0.195498523 -0.813667830 -1.505844690 then, I want to calculate newozone - s(ratio,bs="cr") we define this term (newozone - s(ratio,bs="cr")) is X can i use the coefficient to calculate? -- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3900848.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at 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.
I'm confused... I type .. predict.gam(ozonea,type=s(ratio,bs="cr")) pressure maxtemp s(avetemp) s(ratio) 1 -0.0459102290 -0.185178463 0.263358446 -0.164558673 2 -0.0286464652 -0.194731320 0.199315027 0.727823293 3 0.0478073459 -0.013227033 0.002228896 0.342373202 4 -0.0089164494 0.082301539 -0.037331159 -0.067260889 5 0.0675373617 0.024984396 -0.047067786 -0.357569389 6 0.0823348735 0.101407254 -0.075884852 -1.485036738 7 -0.0977015204 0.177830112 -0.094755158 0.236575309 8 -0.0903027645 0.225594398 -0.113346667 0.435141242 9 0.0206785742 0.187382969 -0.066346157 -0.256133513 10 -0.1371615520 0.101407254 -0.131656887 0.145057584 11 -0.0477674066 -0.181001505 0.260279546 0.180513043 12 -0.0921599421 -0.009050075 -0.020511366 0.281470433 13 0.0681464361 -0.219212934 0.335348247 0.270813178 ...... I want to show s(ratio,bs="cr") term, and show the warning message Warning message: In predict.gam(ozonea, type = s(ratio, bs = "cr")) : Unknown type, reset to terms. I don't understand what does it mean. By the way, i use log-link function, should I convert log-link function of newozone to fitted value of newozone? 1. log(newozone) - s(ratio,bs="cr") = x then exp(x) 2. exp(newozone) - s(ratio,bs="cr") =X then x which one is correct? so confused -- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3901842.html Sent from the R help mailing list archive at Nabble.com.
On Thu, 2011-10-13 at 08:20 -0700, pigpigmeow wrote:
I'm confused... I type .. predict.gam(ozonea,type=s(ratio,bs="cr"))
That is not a valid 'type'; normally you'd use `type = "terms"` or `type = "iterms"`, depending on whether you want (any) standard errors to include the uncertainty about the intercept.
pressure maxtemp s(avetemp) s(ratio) 1 -0.0459102290 -0.185178463 0.263358446 -0.164558673 2 -0.0286464652 -0.194731320 0.199315027 0.727823293 3 0.0478073459 -0.013227033 0.002228896 0.342373202 4 -0.0089164494 0.082301539 -0.037331159 -0.067260889 5 0.0675373617 0.024984396 -0.047067786 -0.357569389 6 0.0823348735 0.101407254 -0.075884852 -1.485036738 7 -0.0977015204 0.177830112 -0.094755158 0.236575309 8 -0.0903027645 0.225594398 -0.113346667 0.435141242 9 0.0206785742 0.187382969 -0.066346157 -0.256133513 10 -0.1371615520 0.101407254 -0.131656887 0.145057584 11 -0.0477674066 -0.181001505 0.260279546 0.180513043 12 -0.0921599421 -0.009050075 -0.020511366 0.281470433 13 0.0681464361 -0.219212934 0.335348247 0.270813178 ...... I want to show s(ratio,bs="cr") term, and show the warning message Warning message: In predict.gam(ozonea, type = s(ratio, bs = "cr")) : Unknown type, reset to terms. I don't understand what does it mean.
So if the above comments I made are not sufficient to explain this, read ?predict.gam and see what the allowed options for the 'type' argument are. G
By the way, i use log-link function, should I convert log-link function of newozone to fitted value of newozone? 1. log(newozone) - s(ratio,bs="cr") = x then exp(x) 2. exp(newozone) - s(ratio,bs="cr") =X then x which one is correct? so confused -- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3901842.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at 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/ UK. WC1E 6BT. [w] http://www.freshwaters.org.uk %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
On Oct 13, 2011, at 11:20 AM, pigpigmeow wrote:
I'm confused... I type .. predict.gam(ozonea,type=s(ratio,bs="cr")) pressure maxtemp s(avetemp) s(ratio) 1 -0.0459102290 -0.185178463 0.263358446 -0.164558673 2 -0.0286464652 -0.194731320 0.199315027 0.727823293 3 0.0478073459 -0.013227033 0.002228896 0.342373202 4 -0.0089164494 0.082301539 -0.037331159 -0.067260889 5 0.0675373617 0.024984396 -0.047067786 -0.357569389 6 0.0823348735 0.101407254 -0.075884852 -1.485036738 7 -0.0977015204 0.177830112 -0.094755158 0.236575309 8 -0.0903027645 0.225594398 -0.113346667 0.435141242 9 0.0206785742 0.187382969 -0.066346157 -0.256133513 10 -0.1371615520 0.101407254 -0.131656887 0.145057584 11 -0.0477674066 -0.181001505 0.260279546 0.180513043 12 -0.0921599421 -0.009050075 -0.020511366 0.281470433 13 0.0681464361 -0.219212934 0.335348247 0.270813178 ...... I want to show s(ratio,bs="cr") term, and show the warning message Warning message: In predict.gam(ozonea, type = s(ratio, bs = "cr")) : Unknown type, reset to terms. I don't understand what does it mean. By the way, i use log-link function, should I convert log-link function of newozone to fitted value of newozone? 1. log(newozone) - s(ratio,bs="cr") = x then exp(x)
Perhaps. I do not see 'newozone' defined anywhere above and can only speculateo speculate. I'm pretty sure Wood's uncopied comment was intended to remind you that differences in the log (regression) scale are equivalent to ratios on the measured scales. So either : exp( diff(coef, log(measured_or_external_baseline) ) ) # OR exp(coef)/measured_or_external_baseline How that is properly calculated in pkg:mgcv is a separate issue and Simpson has suggested you read the help pages more carefully.
2. exp(newozone) - s(ratio,bs="cr") =X then x
Probably incorrect.
which one is correct? so confused
David Winsemius, MD West Hartford, CT
pre<-predict(ozonea,groupA,type="terms",terms=NULL,newdata.guaranteed=FALSE,na.action=na.pass) yeah! but I don't know how to only show the value of s(ratio,bs="cr") -- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3903753.html Sent from the R help mailing list archive at Nabble.com.
predict.gam is returning a matrix with a named column for each term. Select the appropriate column. Example below library(mgcv) n<-200;sig <- 2 dat <- gamSim(1,n=n,scale=sig) b <- gam(y~s(x0)+s(I(x1^2))+s(x2)+offset(x3),data=dat) newd <- data.frame(x0=(0:30)/30,x1=(0:30)/30,x2=(0:30)/30,x3=(0:30)/30) pred <- predict.gam(b,newd,type="terms") pred[,3] ## extracts s(x2) pred[,"s(x2)"] ## same
On 14/10/11 04:54, pigpigmeow wrote:
pre<-predict(ozonea,groupA,type="terms",terms=NULL,newdata.guaranteed=FALSE,na.action=na.pass) yeah! but I don't know how to only show the value of s(ratio,bs="cr") -- View this message in context: http://r.789695.n4.nabble.com/Question-about-GAMs-tp3900848p3903753.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at 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.
Simon Wood, Mathematical Science, University of Bath BA2 7AY UK +44 (0)1225 386603 http://people.bath.ac.uk/sw283