Conditional densities, on one single graph

With Stéphane Tufféry we’ve been working on credit scoring1 and we’ve been using the popular german credit dataset,

> myVariableNames <- c("checking_status","duration","credit_history",
+ "purpose","credit_amount","savings","employment","installment_rate",
+ "personal_status","other_parties","residence_since","property_magnitude",
+ "age","other_payment_plans","housing","existing_credits","job",
+ "num_dependents","telephone","foreign_worker","class")
> credit = read.table(
+ "",
+ header=FALSE,col.names=myVariableNames)
> credit$class <- credit$class-1

We wanted to get a nice code to produce a graph like the one below,

Yesterday, Stéphane came up with the following code, that can easily be adapted

> library(RColorBrewer)
> CL=brewer.pal(6, "RdBu")
> varQuanti = function(base,y,x)
+ {
+ layout(matrix(c(1, 2), 2, 1, byrow = TRUE),heights=c(3, 1))
+	par(mar = c(2, 4, 2, 1))
+	base0 <- base[base[,y]==0,]
+	base1 <- base[base[,y]==1,]
+	xlim1 <- range(c(base0[,x],base1[,x]))
+	ylim1 <- c(0,max(max(density(base0[,x])$y),max(density(base1[,x])$y)))
+	plot(density(base0[,x]),main=" ",col=CL[1],ylab=paste("Density of ",x),
+		 xlim = xlim1, ylim = ylim1 ,lwd=2)
+	par(new = TRUE)
+	plot(density(base1[,x]),col=CL[6],lty=1,lwd=2,
+		 xlim = xlim1, ylim = ylim1,xlab = '', ylab = '',main=' ')
+	legend("topright",c(paste(y," = 0"),paste(y," = 1")),
+		   lty=1,col=CL[c(1,6)],lwd=2)
+	texte <- c("Kruskal-Wallis'Chi² = \n\n",
+       round(kruskal.test(base[,x]~base[,y])$statistic*1000)/1000)
+	text(xlim1[2]*0.8, ylim1[2]*0.5, texte,cex=0.75)
+	boxplot(base[,x]~base[,y],horizontal = TRUE,xlab= y,col=CL[c(2,5)])
> varQuanti(credit,"class","duration")

The code is not complex, but since I usually waste a lot of time on my graphs, I will try to upload more frequently short posts, dedicated to graphs, in R (without ggplot).

1.for a chapter on statistical learning in the forthcoming Computational Actuarial Science with R

2 thoughts on “Conditional densities, on one single graph”

  1. The conditional densities and the boxplots show the conditional distribtuion of duration given class, i.e., P(duration | class). However, in this situation one may argue that P(class | duration) would actually be more interesting. Two types of displays that are easily available in R to show this are spinograms (function spineplot also available through the formula plot method) and conditional density plots (function cdplot). Both want a factor response, hence I added the following to your data processing:

    credit$class <- factor(credit$class, levels = 0:1, labels = c("good", "bad"))

    The spineplot() uses hist() to cut the data into intervals and then plots the conditional proportion of P(class | duration) against P(duration), i.e., intervals with more observations are wider than intervals with fewer observations. This is the same idea as in mosaic plots.

    plot(class ~ duration, data = credit)

    Or using selected duration intervals and switching the order of response levels:

    plot(class ~ duration, data = credit, ylevels = c("bad", "good"),
    breaks = c(4, 6, 11, 12, 18, 24, 36, 72))

    This shows directly that among the very short credits (less than half a year) only 10% were bad while among the long credits (more than three years) around 50% were bad.

    Alternatively, the cdplot() obtains smoothed density() functions for P(duration | class) and P(duration) and then computes P(class | duration) = P(duration | class) * P(class) / P(duration):

    cdplot(class ~ duration, data = credit, ylevels = c("bad", "good"))

    This shows duration directly on the x-axis, rather than P(duration). This can be both an advantage and a disadvantage. As there are only a handful observations beyond 4 years of duration, the conditional density curves are somewhat eratic. It would be safer to just look at

    cdplot(class ~ duration, data = credit, ylevels = c("bad", "good"),
    xlim = c(4, 48))

    This leads to much the same conclusions as the spinogram above. Personally, I prefer the spinogram because it automatically "downweights" areas of the x-axis with few observations.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.