CODEWITHSUNDEEP
rggplot2
Zach
Zach

Reputation: 996

loess and glm plotting with ggplot2

I am trying to plot the model predictions from a binary choice glm against the empirical probability using data from the titanic. To show differences across class and sex I am using faceting, but I have two things things I can't quite figure out. The first is that I'd like to restrict the loess curve to be between 0 and 1, but if I add the option ylim(c(0,1)) to the end of the plot, the ribbon around the loess curve gets cut off if one side of it is outside the bound. The second thing I'd like to do is draw a line from the minimum x-value (predicted probability from the glm) for each facet, to the maximum x-value (within the same facet) and y = 1 so as to show glm predicted probability.

loess-titanic

#info on this data http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3info.txt
load(url('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.sav'))
titanic <- titanic3[ ,-c(3,8:14)]; rm(titanic3)
titanic <- na.omit(titanic) #probably missing completely at random
titanic$age <- as.numeric(titanic$age)
titanic$sibsp <- as.integer(titanic$sibsp)
titanic$survived <- as.integer(titanic$survived)

training.df <- titanic[sample(nrow(titanic), nrow(titanic) / 2), ]
validation.df <- titanic[!(row.names(titanic) %in% row.names(training.df)), ]


glm.fit <- glm(survived ~ sex + sibsp + age + I(age^2) + factor(pclass) + sibsp:sex,
               family = binomial(link = "probit"), data = training.df)

glm.predict <- predict(glm.fit, newdata = validation.df, se.fit = TRUE, type = "response")

plot.data <- data.frame(mean = glm.predict$fit, response = validation.df$survived,
                        class = validation.df$pclass, sex = validation.df$sex)

require(ggplot2)
ggplot(data = plot.data, aes(x = as.numeric(mean), y = as.integer(response))) + geom_point() +
       stat_smooth(method = "loess", formula = y ~ x) +
       facet_wrap( ~ class + sex, scale = "free") + ylim(c(0,1)) + 
       xlab("Predicted Probability of Survival") + ylab("Empirical Survival Rate")

Upvotes: 4

Views: 4581

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226247

The answer to your first question is to use coord_cartesian(ylim=c(0,1)) instead of ylim(0,1); this is a moderately FAQ.

For your second question, there may be a way to do it within ggplot but it was easier for me to summarize the data externally:

g0 <- ggplot(data = plot.data, aes(x = mean, y = response)) + geom_point() +
             stat_smooth(method = "loess") +
             facet_wrap( ~ class + sex, scale = "free") + 
             coord_cartesian(ylim=c(0,1))+
             labs(x="Predicted Probability of Survival",
                  y="Empirical Survival Rate")

(I shortened your code slightly by eliminating some default values and using labs.)

ss <- ddply(plot.data,c("class","sex"),summarise,minx=min(mean),maxx=max(mean))
g0 + geom_segment(data=ss,aes(x=minx,y=minx,xend=maxx,yend=maxx),
                  colour="red",alpha=0.5)

Upvotes: 2

Related Questions