[R] mgcv

Gavin Simpson gavin.simpson at ucl.ac.uk
Thu Jun 10 13:29:23 CEST 2010


On Mon, 2010-06-07 at 10:25 -0700, Dipa Hari wrote:
> 
> Hello Sir,
> I am using mgcv package for my data. 
> My model isy~x1+f(x2),I want to find out the function f(x2). 
> Following is the code.
> 
> sm1=gam(y~x1+s(x2),family=binomial, f) 
> summary(sm1) 
> plot(sm1,residuals=TRUE, xlab="AGE",pch=20) 
>  
> In this plot I am getting S(x2,1.93) on y axixs
> How should I get the function for x2 from this plot.or Is there
> anyother procedureinR to get this function or value for that
> particular function.for example f(x2)= log(x2) so from that plot how
> can I get this kind of formula for x2 variable?

## evaluate smooth at 200 equally spaced points across range of x2
n <- 200
ndat <- with(f, data.frame("x2" = seq(min(x2), max(x2), length = n))
pred <- predict(sm1, newdata = ndat, type = "terms", terms = "s(x2)")

## or for the observed data
pred <- predict(sm1, type = "terms", terms = "s(x2)")

> How can I find variance covariance matrix for that model

That can be done with the standard vcov() extractor method for "gam"
objects:

vcov(sm1)

> and confidence interval for that model ?

Model or smooth s(x2)? If the latter, add se.fit = TRUE to the predict
calls above. If the former;

pred <- predict(sm1, se.fit = TRUE)

or for predictions on scale of your response:

pred <- predict(sm1, type = "response", se.fit = TRUE)

And you can add newdata = ndat to either call if you want predictions
for the newdata points we generated earlier.

>  I tried delta method but I am not getting result.
> Hope you understand my question? 
> Thanks

HTH

G

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 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
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%



More information about the R-help mailing list