Thanks for the very helpful reply. Some comments inline.
On 9/18/2013 8:53 PM, Dennis Murphy wrote:
> Hi Michael:
>
> Some questions:
>
> - Is it possible, and if so, how, to plot the same data and fitted smooths
> on the logit
> scale, i.e., the linear predictor for the binomial glm?
> Yes, but not through stat/geom_smooth directly - the geom only
> provides a simple default mechanism. You can create a data frame using
> predict.glm() with the default type, set se.fit = TRUE and then use
> geom_line() and geom_ribbon() in ggplot2. See below.
Hmm. I thought that ggplot had the facility to apply a transformation,
and found coord_trans,
which I think does what I want, more or less (except that geom_point
doesn't work).
logit <- function(x) log(x)/log(1-x)
ggplot(Titanicp, aes(age, as.numeric(survived)-1, color=sex)) +
stat_smooth(method="glm", family=binomial, formula=y~x,
alpha=0.2, size=2, aes(fill=sex)) +
# geom_point(position=position_jitter(height=0.02, width=0), size=1.5) +
coord_trans(y="logit") + labs(x = "Age", y = "Estimated logit")
> . We first need to get the predicted logits with corresponding SE
> estimates into a data frame along with age and sex, and we should be
> ready to go:
...
With this setup, I'd be happy to plot the observations on the logit
scale jittered around
some reasonable values, say \pm 1.5. However my attempt to do this
doesn't work
and I'm not sure why.
mod <- glm(survived ~ age*sex, family=binomial, data=Titanicp)
modp <- cbind(Titanicp[, c("survived", "sex", "age")],
predict(mod, se.fit = TRUE))
modp$obs <- c(-1.5, 1.5)[modp$survived]
# Plot predicted logits with corresponding Wald CIs
p <- ggplot(modp, aes(x = age, y = fit, color = sex)) +
geom_line(size = 2) +
geom_ribbon(aes(ymin = fit - 1.96 * se.fit,
ymax = fit + 1.96 * se.fit,
fill = sex), alpha = 0.2,
color = "transparent") +
labs(x = "Age", y = "Estimated logit")
p + geom_point(y=modp$obs, position=position_jitter(height=0.02, width=0))
>p + geom_point(y=modp$obs, position=position_jitter(height=0.02, width=0))
Error: position_jitter requires the following missing aesthetics: y
>p + geom_point(y=modp$obs, position=position_jitter(y=modp$obs, height=0.02,
>width=0))
Error in position_jitter(y = modp$obs, height = 0.02, width = 0) :
unused argument (y = modp$obs)
>
> Dennis
>
>> i.e., glm() is quite happy to fit the model survived ~ age+sex
>> in the binomial family, and gives the same predicted probabilities and
>> logits.
>>
>> install.packages("vcdExtra")# data from the most recent version,
>> vcdExtra_0.5-11
>> data(Titanicp, package="vcdExtra")
>> str(Titanicp)
>>
>> 'data.frame': 1309 obs. of 6 variables:
>> $ pclass : Factor w/ 3 levels "1st","2nd","3rd": 1 1 1 1 1 1 1 1 1 1 ...
>> $ survived: Factor w/ 2 levels "died","survived": 2 2 1 1 1 2 2 1 2 1 ...
>> $ sex : Factor w/ 2 levels "female","male": 1 2 1 2 1 2 1 2 1 2 ...
>> $ age : num 29 0.917 2 30 25 ...
>> $ sibsp : num 0 1 1 1 1 0 1 0 2 0 ...
>> $ parch : num 0 2 2 2 2 0 0 0 0 0 ...
>>
>> require(ggplot2)
>> # remove missings on age
>> Titanicp <- Titanicp[!is.na(Titanicp$age),]
>>
>> ggplot(Titanicp, aes(age, as.numeric(survived)-1, color=sex)) +
>> stat_smooth(method="glm", family=binomial, formula=y~x, alpha=0.2,
>> size=2, aes(fill=sex)) +
>> geom_point(position=position_jitter(height=0.02, width=0), size=1.5)
>>
>> # equivalent logistic regression model, survived as a factor
>> mod <- glm(survived ~ age+sex, family=binomial, data=Titanicp)
>> summary(mod)
>>
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept. & Chair, Quantitative Methods
York University Voice: 416 736-2100 x66249 Fax: 416 736-5814
4700 Keele Street Web: http://www.datavis.ca
Toronto, ONT M3J 1P3 CANADA
[[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.
