[R] gam - Y axis probability scale with confidence/error lines

[Patrick Breheny]

Actually, I responded a bit too quickly last time, without really 
reading through your example carefully.  You're fitting a logistic 
regression model and plotting the results on the probability scale.  The 
better way to do what you propose is to obtain the confidence interval 
on the scale of the linear predictor and then transform to the 
probability scale, as in:

x <- seq(0,1,by=.01)
y <- rbinom(length(x),size=1,p=x)
require(gam)
fit <- gam(y~s(x),family=binomial)
pred <- predict(fit,se.fit=TRUE)
yy <- binomial()$linkinv(pred$fit)
l <- binomial()$linkinv(pred$fit-1.96*pred$se.fit)
u <- binomial()$linkinv(pred$fit+1.96*pred$se.fit)
plot(x,yy,type="l")
lines(x,l,lty=2)
lines(x,u,lty=2)

-- 
Patrick Breheny
Assistant Professor
Department of Biostatistics
Department of Statistics
University of Kentucky




On 03/14/2012 01:49 PM, Ben quant wrote:
> That was embarrassingly easy. Thanks again Patrick! Just correcting a
> little typo to his reply. this is probably what he meant:
>
> pred = predict(fit,data.frame(x=xx),type="response",se.fit=TRUE)
> upper = pred$fit + 1.96 * pred$se.fit
> lower = pred$fit - 1.96 * pred$se.fit
>
> # For people who are interested this is how you plot it line by line:
>
> plot(xx,pred$fit,type="l",xlab=fd$getFactorName(),ylab=ylab,ylim=
> c(min(down),max(up)))
> lines(xx,upper,type="l",lty='dashed')
> lines(xx,lower,type="l",lty='dashed')
>
> In my opinion this is only important if the desired y axis is different
> than what plot(fit) gives you for a gam fit (i.e fit <-
> gam(...stuff...)) and you want to plot the confidence intervals.
>
> thanks again!
>
> Ben
>
> On Wed, Mar 14, 2012 at 10:39 AM, Patrick Breheny
> <patrick.breheny at uky.edu <mailto:patrick.breheny at uky.edu>> wrote:
>
>     The predict() function has an option 'se.fit' that returns what you
>     are asking for.  If you set this equal to TRUE in your code:
>
>     pred <- predict(fit,data.frame(x=xx),__type="response",se.fit=TRUE)
>
>     will return a list with two elements, 'fit' and 'se.fit'.  The
>     pointwise confidence intervals will then be
>
>     pred$fit + 1.96*se.fit
>     pred$fit - 1.96*se.fit
>
>     for 95% confidence intervals (replace 1.96 with the appropriate
>     quantile of the normal distribution for other confidence levels).
>
>     You can then do whatever "stuff" you want to do with them, including
>     plot them.
>
>     --Patrick
>
>
>     On 03/14/2012 10:48 AM, Ben quant wrote:
>
>         Hello,
>
>         How do I plot a gam fit object on probability (Y axis) vs raw
>         values (X
>         axis) axis and include the confidence plot lines?
>
>         Details...
>
>         I'm using the gam function like this:
>         l_yx[,2] = log(l_yx[,2] + .0004)
>         fit<- gam(y~s(x),data=as.data.frame(__l_yx),family=binomial)
>
>         And I want to plot it so that probability is on the Y axis and
>         values are
>         on the X axis (i.e. I don't want log likelihood on the Y axis or
>         the log of
>         my values on my X axis):
>
>         xx<- seq(min(l_yx[,2]),max(l_yx[,2]__),len=101)
>         plot(xx,predict(fit,data.__frame(x=xx),type="response"),__type="l",xaxt="n",xlab="Churn"__,ylab="P(Top
>         Performer)")
>         at<- c(.001,.01,.1,1,10)  #<-------------- I'd also like to
>         generalize
>         this rather than hard code the numbers
>         axis(1,at=log(at+ .0004),label=at)
>
>         So far, using the code above, everything looks the way I want.
>         But that
>         does not give me anything information on
>         variability/confidence/__certainty.
>         How do I get the dash plots from this:
>         plot(fit)
>         ...on the same scales as above?
>
>         Related question: how do get the dashed values out of the fit
>         object so I
>         can do 'stuff' with it?
>
>         Thanks,
>
>         Ben
>
>         PS - thank you Patrick for your help previously.
>
>                 [[alternative HTML version deleted]]
>
>         ________________________________________________
>         R-help at r-project.org <mailto:R-help at r-project.org> mailing list
>         https://stat.ethz.ch/mailman/__listinfo/r-help
>         <https://stat.ethz.ch/mailman/listinfo/r-help>
>         PLEASE do read the posting guide
>         http://www.R-project.org/__posting-guide.html
>         <http://www.R-project.org/posting-guide.html>
>         and provide commented, minimal, self-contained, reproducible code.
>
>