3
$\begingroup$

The Hosmer-Lemeshow test has some inconveniences and I am well aware of it. But assuming that I want to apply it then I read that (e.g. Applied logistic regression, by Hosmer and Lemeshow) the test statistic has $g-2$ degrees if freedom ($g$ is the number of partitions used to define the test statistic).

Some of my colleagues however argue that this number of degrees of freedom may be dependent on whether you compute the test statistic on the training sample or on a validation sample. According to them the degrees of freedom to be used on a validation set would be $g$ in stead of $g-2$ because on the validation set you do not estimate any parameters.

Can anyone explain how I should determine the degrees of freedom for (1) the training set and (2) The validation set.

EDIT 8/10/2016

Because I found the existing answers ''confusing'' I added my own answer, I used simulations to make my point clear.

EDIT,

referring to the answer of @jwimberley and the comments below it:

in their 1980 paper, Hosmer and Lemeshow have proven that the test statistic (for a partition with $g$ groups) is $\chi^2(g-p-1) + \sum_{i=1}^p \lambda_i \chi^2(1)$.

One can read that paper and you will find out that the term $\sum_{i=1}^p \lambda_i \chi^2(1)$ is a consequence of the fact that the partition is defined on predicted probabilities and therefore ''random''. The first term can be explained by a similar reasoning as Pearson's GOF test.

In a second step (see 1980 paper) Hosmer and Lemeshow show by simulations that the term $\sum_{i=1}^p \lambda_i \chi^2(1)$ can be approximated by a $\chi^2(p-1)$, and combining this we find that $\chi^2(g-p-1) + \sum_{i=1}^p \lambda_i \chi^2(1)$ is a $\chi^2$ with $g-p-1+p-1=g-2$ degrees of freedom.

All these things will be confirmed by people that read that paper.

The $g-2$ df of the HL statistic are widely known, so any simulation should be able to reproduce that. If not then there is either a problem with the simulation or with the HL paper. I have analysed the HL paper and I think there things are based on mathematical theorems and sound simulations.

I would like to find out how an out-of-sample test (or a test on a validation set) would change the results of Hosmer and Lemeshow, i.e. where in their proof/simulation would there be a difference using a validation set ?

EDIT 30/9/2016

@jwimberley

If you do the Hosmer-Lemeshow test on a validation sample, then would I could expect is that in $\chi^2(g-p-1) + \sum_{i=1}^p \lambda_i \chi^2(1)$ the $p$ in the first term $\chi^2(g-p-1)$, where the $p$ is the consequence of the estimation of $p$ parameters, well for a validation sample there may be an difference on that term because in that validation sample you do not estimate $p$ parameters.

However, for the second term $\sum_{i=1}^p \lambda_i \chi^2(1)$, a term that is there because of the use of predicted (random) probabilities to partition your validation set, well that term should also be there for the validation set. This is what I meant when I said that the degrees of freedom where unusual in this thread: Dividing a sample based on the value of y would be problematic?.

In your first comment and other comments below the answer in this linked thread, you denied that as you can read there.

EDIT 7/10/2016, @jwimberley enter image description here

$\endgroup$
  • $\begingroup$ Have you seen stats.stackexchange.com/questions/167483/… $\endgroup$ – Christoph Hanck Sep 29 '16 at 7:56
  • $\begingroup$ @Christoph Hanck: did you check who posted the question that you refer to? :-) $\endgroup$ – user83346 Sep 29 '16 at 8:03
  • $\begingroup$ No, sorry :-). I just quickly remembered there was a popular question about that topic, that I myself know little about. $\endgroup$ – Christoph Hanck Sep 29 '16 at 8:04
1
$\begingroup$

Note: I have no readily available source for all of the following, but I have simulations demonstrating everything.

It is a general rule that fitting a regression overall always overfits the training by some amount, and that various goodness-of-fit tests/statistics don't behave the same way in the training sample as in the validation sample. Also, the null hypothesis $H_0$ when applying a goodness-of-fit test is different between the two: in the goodness-of-fit test on the training sample, $H_0$ is generally something like "the regression model (e.g. y~x, y~x+I(x^2), etc.) is appropriate." It doesn't require that the regression parameters themselves are exact. On the other hand, the null hypothesis when applying the test to a validation sample is generally that the parameters are exactly correct, since they haven't been tuned to the sample at hand.

A familiar case in point is linear regression. The mean-square-error (MSE) is always lower in the training sample than a validation sample, for the obvious reason that it has been explicitly minimized in the training sample. Consider the following simulation:

n <- 10
a <- 0.1
b <- 0.5
sigma <- 0.1
nums <- 10000
mse <- data.frame(fit=numeric(nums),raw=numeric(nums))
fitmses <- c()
for (i in 1:nums) {
  x <- rnorm(n)
  pred <- a+b*x
  y <- rnorm(n,pred,sigma)
  mse$raw[i] <- mean((y-pred)^2)
  fit <- lm(y~x)
  pred <- predict(fit)
  mse$fit[i] <- mean((y-pred)^2)
}
library(reshape2)
library(ggplot2)
mmse <- melt(mse)
ggplot(mmse,aes(x=value))+
  facet_wrap(~variable,scales = "free_x")+
  geom_histogram()+
  theme_bw()

In this simulation, the true relationship between y and x is $y = 0.1 + 0.5x + \epsilon$ where $\epsilon \sim N(0,0.1)$. Ten thousand simulations with randomly generated x and y are done comparing two quantities: the MSE of fits performed to the sample and the MSE using the correct formula. The two quantities appear somewhat similar:

enter image description here

However, the MSE computed with the true coefficient values is always higher than the the fit value (as it should be). You can easily work out that on the validation sample the MSE uses the exact values, times 10 (the number $n$ of entries in the sample) and divided by $\sigma^2$, follows a chi-square distribution with $n=10$ degrees of freedom:

$$ \sum_i \frac{(\hat y_i - y_i)^2}{\sigma^2} = \sum_i\left(\frac{\epsilon}{\sigma} \right)^2 $$

where $\epsilon/\sigma$ is standard normal, and each term is independent (since a fit hasn't been performed to this dataset). Directly showing this,

h <- hist(mse$raw*n/sigma^2,50)
x <- seq(0,30,0.1)
y <- dchisq(x,df=n)
scale <- (h$count/h$density)[1]
lines(x,y*scale)

enter image description here

Maybe you can do me one better and show this using fitdistr. Remember that this distribution is demonstrating the null hypothesis: that when the regression parameters are exactly correct, this quantity will follow a chi-square distribution with $n$ degrees of freedom. In practice, if a and b have been determined from a fit to a separate training sample, they won't be exactly correct. Then the null hypothesis is invalid. The effective number of degrees of freedom will actually be higher than $n$, leading to larger values of the MSE and low $p$-values from the chi-square test that could reject the null hypothesis.

The MSE from the fitted values is lower, since it has been explicitly minimized, and the rescaled value follows a chi-square distribution with $n-2$ degrees of freedom:

h <- hist(mse$fit*n/sigma^2,50)
x <- seq(0,30,0.1)
y <- dchisq(x,df=n-2)
scale <- (h$count/h$density)[1]
lines(x,y*scale)

enter image description here

Again, if you doubt this you can check it with fitdistr. [Note: I used the known value of sigma in the above in calculating the above distribution. In reality this quantity would be unknown. The lm fit result contains the estimated sigma, which could be used in place, but I won't be redoing the simulation with this]. The values of a and b don't need to be exact, here, because the null hypothesis is that the model $y=a+bx+\epsilon$ is a good representation of the data, not that a and b are exactly correct. The reduction by $2$ in the degrees of freedom related to the two fit parameters. You can try quadratic, cubic, etc. models with $p$ parameters; the reduction in the d.o.f. will increase with $p$, though it won't be exactly $p$ unless $n \gg p$.

In summary: the very act of performing a regression on a training dataset tries to minimize some quantity measuring the discrepancy between predictions and actual values. This means that measures of discrepancy like goodness-of-fit tests have values that are smaller in training sets than in validation sets.

The case of the Hosmer-Lemeshow test for logistic regression is very similar. Here are two small simulations analogous to the above linear regression simulation. In the first, there is an independent variables $x$ and a dependent binary variable $y$. Previous research has shown that the probability that $y$ is true is $\pi = 1/(1+\exp(-[0.1+0.5x]))$. Assuming this is exactly true (the null hypothesis), the following code simulates 10 thousand cross-checks where this prediction is applied to new data and checked against the known values $y$, at each step calculating the Hosmer-Lemeshow statistic. The result follows a chi-square distribution with 10 degrees of freedom -- one for each partition used in the test:

library(ResourceSelection)
n <- 500; chisqs <- c()
for (i in 1:10000) {
    x <- rnorm(n)
    p <- plogis(0.1+0.5*x)
    y <- rbinom(n,1,p)
    x2 <- hoslem.test(y,p)$statistic
    chisqs <- cbind(chisqs,x2);
}
h <- hist(chisqs,50)
xs <- seq(0,25,0.01)
ys <- dchisq(xs,df=10)
lines(xs,ys*h$counts/h$density)

enter image description here

On the other hand, if you use the values of $y$ to make a fit and create the predictions, you'll clearly do better because you've peeked at the answers during the minimization. In this case the degrees of freedom in the distribution of the HL statistic are 8:

library(ResourceSelection)
n <- 500; chisqs <- c()
for (i in 1:10000) {
    x <- rnorm(n)
    p <- plogis(0.1+0.5*x)
    y <- rbinom(n,1,p)
    fit <- glm(y~x,family="binomial")
    x2 <- hoslem.test(y,fitted(fit))$statistic
    chisqs <- cbind(chisqs,x2);
}
h <- hist(chisqs,50)
xs <- seq(0,25,0.01)
ys <- dchisq(xs,df=8)
lines(xs,ys*h$counts/h$density)

enter image description here

This is exactly equivalent to the difference in degrees of freedom for the MSE distributions. Since linear regression is more popular, you're more likely to find sources in the literature verifying what I've claimed and shown here in the linear regression use-case.

In general, you should determine the number of degrees of freedom in the null hypothesis by using simulations like the above. When the number of fit parameters is greater than 2, or when $n$ is not much much greater than $p$, things can start breaking down, and $n-p$ is unlikely to be a good estimate of the degrees of freedom in the training sample. $n$ should always be the number of d.o.f. in the validation sample, however (in the null hypothesis that the regression taken from a different training sample is exactly correct).

EDIT

The number of degrees of freedom is still $g-2$ when doing a logistic regression with a quadratic model with three parameters, contradicting an incorrect claim I made in the comments (and hinted at in this answer where I associated the 2 in $g-2$ with the two fit parameters). In simulation,

library(ResourceSelection)
n <- 500
chisqs <- c()
a <- 0.1
b <- 0.5
c <- 0.2
for (i in 1:10000) {
  x2 <- rnorm(n)
  p2 <- plogis(a+b*x2+c*x2^2)
  y2 <- rbinom(n,1,p2)
  fit <- glm(y2~x2+I(x2^2),family="binomial")
  stat <- hoslem.test(y2,fitted(fit))$statistic
  chisqs <- cbind(chisqs,stat)
}
h <- hist(chisqs,50)
xs <- seq(0,25,0.01)
ys <- dchisq(xs,df=8)
lines(xs,ys*h$counts/h$density) 

enter image description here

So, I apologize for this mistake. It does not take away from the rest of the answer, however.

SECOND EDIT

Here is the same quadratic model, on a validation sample, as requested:

library(ResourceSelection)
n <- 500
chisqs <- c()
a <- 0.1
b <- 0.5
c <- 0.2
for (i in 1:10000) {
  x2 <- rnorm(n)
  p2 <- plogis(a+b*x2+c*x2^2)
  y2 <- rbinom(n,1,p2)
  stat <- hoslem.test(y2,p2)$statistic
  chisqs <- cbind(chisqs,stat)
}
h <- hist(chisqs,50)
xs <- seq(0,25,0.01)
ys <- dchisq(xs,df=10)
lines(xs,ys*h$counts/h$density)

enter image description here

FINAL EDIT

As requested, here is one more simulation, in which the parameters a and b are estimated beforehand from a larger training sample with ten times as many entries. Then applying this fit to the validation samples, the do.f. is still close to 10:

library(ResourceSelection)
n <- 500
chisqs <- c()
a <- 0.1
b <- 0.5
c <- 0.2
x <- rnorm(10*n)
p <- plogis(a+b*x+c*x^2)
y <- rbinom(10*n,1,p)
fit <- glm(y~x+I(x^2),family="binomial")
for (i in 1:10000) {
  x2 <- rnorm(n)
  p2 <- plogis(a+b*x2+c*x2^2)
  y2 <- rbinom(n,1,p2)
  stat <- hoslem.test(y2,predict(fit,data.frame(x=x2),type="response"))$statistic
  chisqs <- cbind(chisqs,stat)
}
h <- hist(chisqs,50)
xs <- seq(0,25,0.01)
ys <- dchisq(xs,df=10)
lines(xs,ys*h$counts/h$density)

enter image description here

If the number of data points in the training sample (above, 5000) was decreased, the measured parameters a, b, and c would differ from their exact values by more and the null hypothesis would be violated. The distribution of the HL test statistic would no longer be chi-square and would skew higher (making it look like it had $>g$ d.o.f. if you superimposed a chi-square distribution).

$\endgroup$
  • $\begingroup$ Can you make simulations for this: estimate a logistic regression with 5 parameters and simulate the HL test statistic 'in sample', compare it to a chi-square with $g-2$ degrees of freedom. Can you do that? Because your answer is not about logistic and therefore not about HL and therefore not related to my question, so please simulate HL for a logistic regression, once in sample and once out of sample. That is what my question is about $\endgroup$ – user83346 Sep 29 '16 at 20:16
  • $\begingroup$ If your simulations do what you pretend they do then the in sample HL for a model where you estimated 5 parameters should resemble chi square with g-2 df, so please show me that with your simulations $\endgroup$ – user83346 Sep 29 '16 at 20:21
  • $\begingroup$ @fcop About half my answer is about logistic regression; I started off with linear regression since it simpler and exhibits similar features. I did have a typo where "family="binomial" was left off in the glm logistic regression fit for the second HL simulation I show, so I'm sorry if that caused confusion (the typo wasn't present when I generated the image below the code, so it has no issues). You should be able to modify the simulations as you request, if you would like. I doubt that with 5 estimated parameters the d.o.f. would be $g-2$; it would be somewhere between this and $g-5$. $\endgroup$ – jwimberley Sep 29 '16 at 20:25
  • 1
    $\begingroup$ The H-L test is no longer recommended. It has been replaced by more powerful and actionable tests. So I'm not sure why the discussion is worth this many words. $\endgroup$ – Frank Harrell Oct 9 '16 at 12:00
  • 1
    $\begingroup$ @FrankHarrell It wasn't and I regret it. I'd hoped to make a broader point also true for better tests. E.g. $\sum_i (Y_i - \pi_i)^2$ leads to the the Hosmer-le Cessie test in fit sample and Spiegelhalter test in a validation sample; former requires reduction by quadratic form involving score vector/information matrix to account for overfitting, while the latter is straightforward. Analogous to RSS in linear regression having $n-p$ d.f. in fit sample but $n$ in validation, and H-L test having $g-2$ in fit sample and $g$ in validation (subject of original question before recent title change). $\endgroup$ – jwimberley Oct 9 '16 at 13:43
1
$\begingroup$

In their 1980 paper Hosmer D.W., S. Lemeshow, ''Goodness of fit tests for the multiple logistic regression model'', Communications in Statistics - Theory and Methods, Volume 9, 1980 - Issue 10, the authors have proven that the test statistic (for a partition with $g$ groups) is $\chi^2(g-p-1) + \sum_{i=1}^p \lambda_i \chi^2(1)$.

In a second step they showed, using simulations that the term $\sum_{i=1}^p \lambda_i \chi^2(1)$ is a approximately $\chi^2(p-1)$ and as the sum of (independent) $\chi^2$ is also $\chi^2$ with degrees of freedom equal to the sum of the individual degrees of freedom

they found that (1) their test statistic is exactly $\chi^2(g-p-1) + \sum_{i=1}^p \lambda_i \chi^2(1)$ and (2) their test statistic is approximately $\chi^2(g-2)$

Simulations under the Hosmer-Lemeshow conditions

As they have shown this formally, any simulation that is executed under their necessary conditions should find (approximately) a $\chi^2(8)$ when $g=10$. This is the case for the simulations below (in order to keep the answer ''readable'', I inserted the ''helper'' functions at the botttom of this answer, these have to be executed first).

The code contains comments to point to the major steps in the simulation. It is mainly a loop that is executed $N$ times, in each loop one has:

  1. Draw a training sample
  2. Estimate the logistic model
  3. Predict the probabilities using the estimated coefficients
  4. Compute the HL X2

The helper functions are defined at the end of my answer.

simulateHosmerLemeshowX2<-function(N=5000, sampleSize, b0=-4) {

  x2HL<-vector(mode="numeric", length=N)

  for (i in 1:N) {

    # draw a training sample
    trainSample<-generateSample(sampleSize, b0)

    # train a logistic model
    logisticModel <- glm(y~x1+x2,family="binomial",data=trainSample$sample)

    #predict the probabilities using the estimates
    predictedPs<-predict(logisticModel, newdata = trainSample$sample, type = "response")

    # compute the Hosmer-Lemeshow statistic
    x2HL[i] <- hoslem.test(trainSample$sample$y,predictedPs)$statistic
  }


  p<-plotSimulatedX2(x2HL, N=N,title = "H0: model predicts probabilities well (Hosmer-Lemeshow)", df.chi=8)
  return(list(graph=p, N=N))    
}

simulateHosmerLemeshowX2(sampleSize=500, b0=-4)

One can execute this simulation, you will find a graph like the one below:

The bars represent the simulated test statistic of Hosmer-Lemeshow, the smooth curve is a simulated $\chi^2(8)$. It seems to conform the findings of HL. Note that the mean of the test statistic is close to 8 and that the mean of a chi-square is the number of degrees of freedom. enter image description here

Simulations for ''out-of-sample'' validation

In this section I simulate the out-of-sample situation. The code and one result are below.

The code contains comments to point to the major steps in the simulation. It is mainly a loop that is executed $N$ times, in each loop one has:

  1. Draw a sample
  2. Split it into a training and a validation sample
  3. Estimate the logistic model on the training sample
  4. Predict the probabilities using the estimated coefficients on the validation sample
  5. Compute the HL X2 on the validation sample

The helper functions are defined at the end of my answer.

The graphs shows that the test statistic, compute on the validation sample, is not $\chi^2(8)$, the mean seems to indicate that one might conclde that it is $\chi^2$ with 13 df, but, as indicated in the graph, then the variance is not compatible with a chi-square ?

simulateHosmerLemeshowOutOfSampleX2B<-function(N=5000, sampleSize, b0) {

  x2HL<-vector(mode="numeric", length=N)


  for (i in 1:N) {

    # generate a sample
    fullSample<-generateSample(2*sampleSize, b0)

    # split the sample in two equal parts; one for training and one for validation
    idx<-sample(x = 1:(2*sampleSize), size = sampleSize)
    trainSample<-fullSample
    trainSample$sample <-fullSample$sample[idx,]
    trainSample$truePs <-fullSample$truePs[idx]

    validationSample<-fullSample
    validationSample$sample <-fullSample$sample[-idx,]
    validationSample$truePs <-fullSample$truePs[-idx]

    # train a logistic model on the training sample
    logisticModel <- glm(y~x1+x2,family="binomial",data=trainSample$sample)

    # use the trained model on the validation sample
    predictedPs<-predict(logisticModel, newdata = validationSample$sample, type = "response")

    # compute HL on the validation sample
    x2HL[i] <- hoslem.test(validationSample$sample$y,predictedPs)$statistic
  }


  p<-plotSimulatedX2(x2HL, N=N,title = "H0: model predicts probabilities well, Out of Sample \n Splitted sample", 
                     df.chi=12)

  return(list(graph=p, N=N)) 
}

simulateHosmerLemeshowOutOfSampleX2B(sampleSize=500, b0=-4, N=5000)

enter image description here

Annex: code for the helper functions.

library(ResourceSelection)
library(ggplot2)


generateSample<-function(sampleSize=100, b0=-4, b1=0.5, b2=3) {

  x1<-rnorm(sampleSize,mean=3)
  x2<-rnorm(sampleSize, mean=1)

  p <- 1/(1+exp(-(b0+b1*x1+b2*x2)))
  y <- rbinom(sampleSize,1,p)

  return(list(sample=data.frame(y=y,x1=x1,x2=x2), truePs=p))
}

plotSimulatedX2<-function(x2HL, N, title="", df.chi=8) {

  df1<-data.frame(x=x2HL, type="X2")
  df2<-data.frame(x=rchisq(n = 3*N, df = df.chi), type=paste("Chi^2(", df.chi,")") )
  df<-rbind(df1, df2)

  meanX2<-round(mean(df1$x),digits=2)
  varX2<-round(var(df1$x),digits=2)

  subTit<-paste("mean(X2)=",meanX2, " \nvar(X2)=", varX2,sep="")

  p<-ggplot() +
    geom_histogram(data = df1, aes(x=x, fill=type, y=..density..),binwidth=2, alpha=.5, position="identity") +
    geom_density(data = df2, aes(x=x, fill=type), alpha=.3)+
    xlim(c(0,50))+
    annotate("text",x=40,y=0.05,label=subTit, colour="red")+
    annotate("rect", xmin = 30, xmax = 50, ymin = 0.042, ymax = 0.058, alpha = .2)+
    ggtitle(label = title)+
    theme(plot.title = element_text(size = rel(1.5)))

  return(p)
}
$\endgroup$
0
$\begingroup$

In the original paper, Hosmer and Lemeshow used 8 df (with 10 decile groups) for estimating the model on development data. For validation data, you would use 9 df (df= # of groups - 1), based on what I've seen in literature.

Edit: I haven't got a reference to support the choice of 9 df.

$\endgroup$
  • $\begingroup$ Can you explain more in detail? Because that is not what I can derive from the paper. $\endgroup$ – user83346 Sep 29 '16 at 15:49
  • $\begingroup$ My apologies, I haven't seen the original and have gleaned this from later sources. You might be interested to see recommendations for the same by Paul Alison (support.sas.com/resources/papers/proceedings14/1485-2014.pdf) and also thestatsgeek.com/2014/02/16/… $\endgroup$ – prince_of_pears Sep 29 '16 at 16:10
  • $\begingroup$ Thanks, I will read that in detail but at first glance none of these two references say anything about degrees of freedom in a validation set or did I miss that? $\endgroup$ – user83346 Sep 29 '16 at 16:34
  • $\begingroup$ These ones don't really about validation, only development. You could look at Harrell's book for the other. $\endgroup$ – prince_of_pears Sep 29 '16 at 17:03
  • $\begingroup$ But in your answer you say 8 for training and 9 for validation, what is that based on? If you can not justify it then I will have to downvote. $\endgroup$ – user83346 Sep 29 '16 at 18:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy