Plotting logit-transformed binomial data fails

[Mike Lawrence]

The code below first fits a smoother to binomial data then transforms
the y-axis to logit scale, which should linearize the smooth.
Unfortunately the plot completely fails. Any ideas why?

data("kyphosis", package="rpart")
kyphosis$y = as.numeric(kyphosis$Kyphosis)-1
ggplot(
data = kyphosis
, mapping = aes(
x = Age
, y = y
)
)+
layer(
stat = 'smooth'
, stat_params = list(
method="glm"
, family="binomial"
)
)+
coord_trans(y='logit')

--
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tr.im/mikes_public_calendar

~ Certainty is folly... I think. ~

[jc]

I wanted to know the answer to this as well. I tried extracting the
results of the smooth (..y.. ?) and then scaling. That doesn't work
either.

[hadley wickham]

> I wanted to know the answer to this as well.  I tried extracting the
> results of the smooth (..y.. ?) and then scaling.  That doesn't work
> either.

This helps to show the problem:

kyphosis$y = as.numeric(kyphosis$Kyphosis)-1
ggplot(kyphosis, aes(Age, y)) +
geom_point() +
stat_smooth(method = "glm", family = "binomial")

There are data points at 0 and 1, which become -Inf and Inf under the
logistic transformation.

Hadley

--
http://had.co.nz/

[jc]

But the documentation states that the transformation happens at the
end with coord_trans. In that case the smoothed values that I'm
actually plotting do not contain 1 or 0. Furthermore, I should be
able to grab the result of the smooth and plot just that (..y.. it
seems from the documentation). In that case there's definitely no 1
or 0 and I should be able to transform it. But, it doesn't work.

[hadley wickham]

> But the documentation states that the transformation happens at the
> end with coord_trans.

Correct.

> In that case the smoothed values that I'm
> actually plotting do not contain 1 or 0.

Also correct.

> Furthermore, I should be
> able to grab the result of the smooth and plot just that (..y.. it
> seems from the documentation).

Unfortunately it's not that simple. This is pretty subtle and took me
a while to get right, but the scales have to be trained with the raw
data as well as the output of the statistics - this is what enables
(among other things) fullrange = T to work for stat_smooth.

Hadley

--
http://had.co.nz/

[jc]

So, the short answer is that it's not possible to do the transform on
a value derived from raw data if the transform cannot be applied to
the raw data as well. So... back to getting predicted values
beforehand and plotting them. Unfortunately I lose the nice 95%CI
ribbon.

[hadley wickham]

> So, the short answer is that it's not possible to do the transform on
> a value derived from raw data if the transform cannot be applied to
> the raw data as well.  So... back to getting predicted values
> beforehand and plotting them.  Unfortunately I lose the nice 95%CI
> ribbon.

It's not hard to calculate the ribbon yourself - see geom_ribbon for
how to plot it.

Hadley

--
http://had.co.nz/

[jc]

It occurs to me that it's still unsatisfactory that I can't even
plot ..y.. explicitly. And, why would that need to be trained on the
raw data? I'm specifying that I only want the stat result. It should
word like the histogram data. There isn't even a raw y there.

BTW, plotting the ..y.. result alone fails. (no y specified error)

[hadley wickham]

> It occurs to me that it's still unsatisfactory that I can't even
> plot ..y.. explicitly.  And, why would that need to be trained on the
> raw data?  I'm specifying that I only want the stat result.  It should
> word like the histogram data.  There isn't even a raw y there.

I don't understand what y you are talking about. Could you please
supply a reproducible example?

Hadley

--
http://had.co.nz/

[jc]

I don't know if I can reproduce something that doesn't work. But, a
simple modification to your example. would be...

kyphosis$y = as.numeric(kyphosis$Kyphosis)-1
ggplot(kyphosis, aes(Age, y)) + geom_point() +

+ stat_smooth(aes(y=..y..), method = "glm", family = "binomial")

Error: stat_smooth requires the following missing aesthetics: y


I got ..y.. because it's an output variable in the stat_smooth
documentation. It's supposed to be y-predicted. So, I just would add
an aes(y = ..y..) with the stat_smooth. Given the error then perhaps
the stat_smooth output is something else? It seems like there should
be some other error that occurs if y really is an output of
stat_smooth.

[hadley wickham]

> I got ..y.. because it's an output variable in the stat_smooth
> documentation.  It's supposed to be y-predicted.  So, I just would add
> an aes(y = ..y..) with the stat_smooth.  Given the error then perhaps
> the stat_smooth output is something else?  It seems like there should
> be some other error that occurs if y really is an output of
> stat_smooth.

It's a little complicated - y is both an input to and an output from

[jc]

OK, well, perhaps I should go on the site and suggest the stat_smooth
documentation be updated to reflect the fact that the y output can't
just be accessed like say, density in histogram. :)