4
$\begingroup$

I am trying to calculate the Hosmer-Lemeshow Goodness-of-Fit (GOF) test, following the steps presented in the book "Logistic Regression. A Self- Learning Text" by David Kleinbaum and M. Klein.

The following function calculates all the steps, gives as outputs the chi-squared, degrees of freedom and the p-value, as well as some interesting plots described in the manual, which look fine to me.

However, I get a completely different result when I run this test in R (holsem.test, from ResourceSelection package).

I was surprised I could not find any code online, so maybe this becomes useful in the future. Neither could I see what I am doing wrong. Maybe R uses a very different method?

If you are also keen to find out and see this work, please feel free to join.

Thank you a lot,

EDIT: one thing that I find intriguing is the fixed lenght or size of the deciles. My data, for example, has 5479 observations, which cannot be nicely divided by 10. This may lead to some discrepancies bewteen methods. However, I believe 5479 is a long enough series and ignoring a maximum of 9 observations should not be of great harm.

 % This function performs the Hoslem-Lemeshow Goodness-of-Fit test
    % 06-January-2017 davidnielsen@id.uff.br
    %
    % Inputs
    %       Q: The number of groups or quantiles (normally Q=10).
    %       y: 1D array of observations. Must be a series of {0,1}. 
    %       fit: 1D array of the model's fitted values.
    % Outputs
    %       chi2: The Hosmer-Lemeshow chi-square statistics.
    %       df: Degrees-of-freedom.
    %       p: p-value

    function [chi2,df,p]=hoslem(fit,y,Q)
    [ordems,I]=sort(fit,'descend'); % Sorts the fitted values
    % Calculates the length of each decile.
    % Note that the division of length/10 may not be exact. Therefore, the
    % considered series may differ from original inputs in up to 10
    % observations, at most. This error shold be to non-significant for very
    % long series.
    qlen=round(length(y)/Q)-1;      
    qlens=zeros(Q+1,1); 
    for i=1:Q+1
        qlens(i,1)=1+(i-1)*qlen;    % Determines the thresholds of deciles   
    end
    qlims=zeros(Q+1,1);             % Positions of thresholds
    for i=1:Q+1
        qlims(i,1)=ordems(qlens(i,1));
    end
    plot(1:Q+1,qlims,'o')
    % Couts Observed Cases por Decile (Ocq)
    ocq=zeros(Q,1); t=1; count=0;
    for i=1:max(qlens)-1
        if count==qlen
            t=t+1; count=0;
        end
        ocq(t,1)=ocq(t,1)+(y(I(i)));
        count=count+1;
    end
     % Couts Non-Observed Cases por Decile (Oncq)
    oncq=zeros(Q,1); t=1; count=0;
    for i=1:max(qlens)-1
        if count==qlen
            t=t+1; count=0;
        end
        if (y(I(i)))==0
        oncq(t,1)=oncq(t,1)+1;
        end
        count=count+1;
    end
    figure(1)
    plot(ocq); hold on; plot(oncq)
    legend('Observed cases','Non-observed cases'); hold off
    % Expected cases (Ecq) e noncases (Encq) per Decile
    ecq=zeros(Q,1);
    t=1; count=0;
    for i=1:max(qlens)-1
        if count==qlen
            t=t+1; count=0;
        end
        ecq(t,1)=ecq(t,1)+ordems(i);
        count=count+1;
    end
    figure(2)
    plot(ecq,ocq,'-o'); hold on;
    xlabel('Expected cases'); ylabel('Oobserved cases');
    encq=qlen-ecq;
    % Calculates HL statistics
    chi2= sum(sum(((ocq-ecq).^2)/ecq)) + sum(sum(((oncq-encq).^2)/encq));
    df=Q-2; 
    p=1-chi2cdf(chi2,Q-2);
    end
$\endgroup$
  • $\begingroup$ Although this question asks about code, the underlying issue is statistical. This can remain open, IMO. $\endgroup$ – gung - Reinstate Monica Jan 6 '17 at 21:24
2
$\begingroup$

When there are multiple responses within each combination of values of your predictor variables, you can compare your predicted probabilities directly to the observed outcomes. (You can see an example in my answer here: Test logistic regression model using residual deviance and degrees of freedom.) The tricky issue is what to do when the observations are not grouped.

Hosmer and Lemeshow (1980) suggested that you sort all the predicted probabilities and bin them somehow. Then you could compare the number of 'successes' in each bin to the expected number (the sum of the predicted probabilities) in that bin with a chi-squared test. This sounds good, as far as it goes, but there is no principled (i.e., correct) way to bin the predicted probabilities. There are lots of ways of trying to bin them, but no way is necessarily right and the test is very sensitive to the binning. This has been noted by no less than Hosmer, Hosmer, Le Cessie and Lemeshow (1997, SIM).

With this understanding, we can infer the reason you get different results from running the Hosmer and Lemeshow test on your same model in MATLAB and R: The binning is being done differently and the result is sensitive to the binning. Moreover, the test (in both MATLAB and R) is unreliable and should be ignored.

For a nice overview of this topic, see: Alison, PD (2014). Measures of Fit for Logistic Regression. SAS Global Forum, paper 1485-2014. (pdf)

$\endgroup$
  • $\begingroup$ Thank you. That is interesting. I read indeed that the binning in quantiles could be done either using number of observations or the fitted values themselves. After exploring a bit more, I see now that a change in the number of quantiles (generally 10) also changes the results. $\endgroup$ – ouranos Jan 7 '17 at 0:42
  • $\begingroup$ You're welcome, @DavidMarcolinoNielsen. $\endgroup$ – gung - Reinstate Monica Jan 7 '17 at 0:59

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.