Pseudo-$R^2$ statistics under complex sampling
Thomas Lumley

Dipolar Dark Matter as an Effective Field Theory
    Lavinia Heisenberg
    December 30, 2023
================================================




Model summaries based on the ratio of fitted and null likelihoods have been proposed for generalised linear models, reducing to the familiar R^2 coefficient of determination in the Gaussian model with identity link.  In this note I show how to define the Cox–Snell and Nagelkerke summaries under arbitrary probability sampling designs, giving a design-consistent estimator of the population model summary.   I also show that for logistic regression models under case–control sampling the usual Cox–Snell and Nagelkerke R^2 are not design-consistent, but are systematically larger than would be obtained with a cross-sectional or cohort sample, even in settings where the weighted and unweighted logistic regression estimators are similar or identical. 

Keywords: likelihood, logistic regression, case-control study, sampling weights




§ BACKGROUND

The coefficient of determination for a linear regression model is the proportional reduction in squared prediction error from using the model prediction instead of the mean. That is, if we have n observations, a vector X of p predictors and a model 

    E[Y|X=x]=μ=α+xβ

we define

    R^2=1-∑_i=1^n (Y_i-μ̂_i)^2/∑_i=1^n (Y_i-Y̅)^2.

<cit.>, in an exercise, proposed a definition for binary Y in terms of the likelihood ratio 

    R^2_CS=1-(L(0)  L(β̂))^2/n

where L(0) is the maximised likelihood  for an intercept-only model 

    E[Y]=α

and L(β̂) is the likelihood maximised over α and β. 

This definition reduces to the coefficient of determination in a linear-Normal model with Y∼ N(μ,σ^2) if the likelihood is also maximised over the variance parameter (so -2log L=2+n+nlog(2πσ̂^2)), though not if L(0)/L(β̂) is the likelihood ratio from a model with fixed σ (as is common for generalised linear models).

<cit.> noted that for Bernoulli data the maximum possible value of the likelihood is 1, and so the maximum possible value of R^2_CS is 1-(L(0))^2/n. He defined a rescaled pseudo-R^2

    R^2=. R^2_CS/ (1- L(0))^2/n.

which attains the value 1 when prediction is perfect.  Nagelkerke's R^2 has become popular for logistic regression. 


As maintainer of the survey package for R<cit.>, I have been asked on more than one occasion how to compute Nagelkerke's R^2 under complex sampling. 
In this note I propose a definition in terms of a superpopulation parameter, and a design-based estimator. I show  the estimator consistently estimates the population value and reduces exactly to Nagelkerke's statistic under simple random sampling.  I also compare the proposed definition to the estimates obtained by naive use of the standard formula in unweighted logistic regression under case–control sampling, the standard analysis method in epidemiology and biostatistics.

Code and data for the examples is on <github.com/tslumley/pseudorsq>, and the estimators are available in the `survey' package starting from version 3.31–7.



§ DEFINITION

Following <cit.>, I will assume a finite population of size N is an independent and identically distributed draw of vectors (X, Y) from some (unknown) probability model with density g(x,y), and a probability sample of size n is taken from with known sampling probabilities π_i for the ith individual.  I will write E_p[·] for expectations with respect to the superpopulation model and E_π[·] for expectations with respect to finite sampling, and write w_i for the sampling weights 1/π_i.

The analysis goal is to fit parametric regression models f(y|x; θ) for the marginal densities of Y given X, where the parameter vector θ includes the regression intercept and slopes (α, β) and possibly other nuisance parameters.  I do not assume this parametric family necessarily contains the true model; the aim is to estimate the `least false' parameter θ^*, the value minimising the Kullback–Leibler divergence between f(y|x; θ) and the true g(·).  Where there is a possibility of confusion, I will refer to f(y|x; θ) as the `working model' and the likelihood based on it as the `working likelihood'.  Let θ̃_N be the `census parameter', the value obtained by fitting the model f(y|x; θ) to the whole population.  By standard results on maximum likelihood estimation, θ̃_Np→θ^*. 

The parametric model is fitted to the complex sample by maximising the weighted loglikelihood 

    ℓ̂(θ) = ∑_i∈sample w_i ℓ_i(θ) = ∑_i∈sample w_i log f(y_i|x_i; θ).

The weighted likelihood is unbiased for the population working loglikelihood: 

    E_π[ℓ̂(θ)]=∑_i=1^N ℓ_i(θ)

and its derivative is unbiased for the population working score function. 
Under mild conditions on the superpopulation distribution g(·) and the sampling design, the resulting estimator θ̂ is √(n)-consistent for θ̃_N, and so also for θ^*.   Because the main interest of the analyst is usually in the regression parameters β I will abuse notation by writing ℓ̂(β) for the value of ℓ̂(θ) with any other parameters estimated by maximum weighted likelihood. 




We can now consider the Cox–Snell summary, R^2_CS. First, let ℓ_1(β) be the working loglikelihood for a single random observation from the superpopulation distribution, and define the superpopulation summary ρ^2_CS by

    log (1- ρ^2_CS) =2E_p[ℓ_1(β^*)-ℓ_1(0)].

The ordinary Cox–Snell R^2_CS for an iid sample is an estimator of ρ^2_CS obtained by replacing the expectation with a population or sample average, and under iid sampling the consistency of R^2_CS for ρ^2_CS follows from the law of large numbers, the consistency of  maximum likelihood, and the smoothness of ℓ(·). 


The finite population has been constructed as an iid sample, so I can define the census (population) parameter R̃^2_CS by

    log (1-R̃^2_CS) =2/N(ℓ(β̃)-ℓ(0))=1/N∑_i=1^N (ℓ_i(β̃)-ℓ_i(0))

and obtain a design-based plug-in estimator

    log (1-R̂^2_CS) =2/N̂(ℓ̂(β̂)-ℓ̂(0))=∑_i∈sample w_i (ℓ_i(β̂)-ℓ_i(0))/∑_i∈sample w_i

That is, under complex sampling, the likelihood should be replaced by the weighted (pseudo)likelihood, and the n in the Cox–Snell and Nagelkerke formulas should be replaced by the sum of the weights.   

The analysis so far has not assumed the model is correctly specified.  If the fitted parametric model family f(·; θ) happened to include the true superpopulation distribution g(·), the right-hand side of equation <ref> would be the mutual information between X and Y, or equivalently, the difference between the marginal and conditional entropy of Y.  This relationship to information theory gives another reason to regard ρ^2_CS as a meaningful model summary.  

The design-based Cox–Snell statistic R̂^2_CS can now be rescaled to an upper limit of 1, obtaining a design-based version of the Nagelkerke pseudo-R^2.

    R̂^2=. R̂^2_CS/ (1- exp(ℓ̂(0)))^2/n.





§ COMPARISONS

In this section I compare the proposed estimator either to a known or simulated population value or to the estimator obtained by ignoring the sampling. 



 §.§ A large multistage survey
 
NHANES, the National Health and Nutrition Examination Survey, is a series of large national surveys of the US civilian population conducted by the National Center for Health Statistics (NCHS). The survey has been run continuously since 1999, with data released for each two-year wave.  Each wave samples about 10,000 individuals, in a complex multistage, multiphase design, which is approximated for the purpose of public-use datasets by a two-stage sample.  In the approximate design, two city or county sampling units are taken from each of about 16 geographical strata, and a total of about 10,000 individual are sampled.  The survey oversamples people under 18 and over 60, and also oversamples racial/ethnic minority groups.  We will use data on blood pressure from the 2003–4 and 2005–6 waves of NHANES, in 18323 individuals for whom both dietary data and blood pressure were available<cit.>.

Isolated Systolic Hypertension (ISH) is the form of hypertension most common in older people. It is defined by an elevated systolic blood pressure (>140 mmHg) with normal or low diastolic blood pressure (<90 mmHg), and indicates stiffening of walls of major arteries<cit.>. Following <cit.> we fit a sequence of logistic regression models. The base model  used age as a linear spline with interior knots at 50 and 65 years. Successive models then added race/ethnicity as five categories, gender, a gender by age interaction, and reported dietary sodium intake.  The R^2 statistics for the models are shown in Table <ref>; for both the Cox–Snell and Nagelkerke statistics the design-based estimator has a slightly lower value than the estimator assuming simple random sampling.





 §.§ Case–control sampling

Logistic regression in case–control designs is probably the most frequently-used example of a generalised linear model under complex sampling.  Case–control sampling is used in the study of rare diseases: the small number of individuals with Y=1 are all sampled, but only a small fraction of those with Y=0. The standard analysis of a case–control sample in epidemiology is by unweighted logistic regression; the bias that would be expected from ignoring the sampling is confined to the intercept and the estimate of β is the semiparametric MLE.  The unweighted analysis is sometimes, but not always, substantially more efficient than the weighted analysis. 



  §.§.§ Heuristics

A heuristic analysis is useful.  Suppose that P(Y=1) is very small in the population, and that π_0≪ 1 is chosen to give m controls per case. The efficiency of this design (with unweighted analysis) relative to using the entire population is m/(m+1), which can be chosen close to 1.   Almost all the (Fisher) information in the population is now present in the sample, so the likelihood ratio  between the null model and the model fitted by maximum likelihood will be approximately the same in the sample as in the population.   The ordinary Cox–Snell R-squared  in the population, R̃^2_CS, satisfies 

    log (1-ρ^2_CS) ≈log (1-R̃^2_CS) =2/N(ℓ(β̃)-ℓ(0))

so in the unweighted model in the sample it will  satisfy

    log (1-R^2_CS,sample) ≈2/n(ℓ(β̃)-ℓ(0))

different by a factor of N/n.  That is, case–control sampling will dramatically deflate log(1-R^2_CS) and so inflate R^2_CS relative to prospective sampling from the same population.  In a less-ideal case–control scenario we would expect the sample likelihood ratio to be smaller than the population likelihood ratio — but if it were smaller by  a factor of n/N, to cancel the bias, a case–control design would have no advantage over simple random sampling. 

The Nagelkerke correction will mitigate the sampling bias, because the rescaling factor is also affected by sampling bias and the biases partly cancel. However, 
 the sampling bias is multiplicative on the scale of the log likelihood ratio and the correction is multiplicative on the scale of R^2, a concave function of the log likelihood ratio. By Jensen's inequality, the ordinary Nagelkerke R^2 in the sample will still be biased upwards relative to the population value except when R^2=0 or R^2=1.  The design-based estimator R̂^2 does not have these biases; it will estimate the same quantity under case–control sampling as under prospective sampling. 
 


  §.§.§ Simulation

Table <ref> demonstrates this bias in a simulation. A population of size 10^5 was simulated with a single N(0,1) predictor, X, and a binary outcome Y satisfying

    logit P(Y=1|X=x)=-6+x.

There were 389 cases in the simulated population, and control samples were drawn with 1, 2,5, 10, or 20 controls per case, giving control sampling fractions from 0.4% to 7.8%. As Table <ref> shows, the design-based estimators accurately reproduced the population R^2 statistics under all sampling fractions. The unweighted estimators displayed a substantial upwards bias that decreased with increasing sampling fraction, as expected from our heuristic argument. 


 
 


  §.§.§ A practical example

To illustrate the practical relevance of this bias, consider the case–control study of oesphageal cancer in men from Ille-et-Vilaine, France, reported by <cit.> and later used by <cit.> a study that is close to the heuristic ideal.  There are two publicly-available data sets: one with discrete variables for age (five groups),  tobacco consumption (four groups), and alcohol consumption (four groups), and the other with continuous variables for age and reported alcohol and tobacco consumption.   The study recruited between four and five controls per case; the control sampling fraction is not given explicitly but can be estimated from population data to be about 1/440. 

For both data sets, a model with main effects of each variable fits reasonably well. We consider this base model and a model adding a linear by linear interaction between  alcohol and tobacco consumption.  When using the grouped data, the estimated variances of β̂ are very similar for the unweighted and design-weighted approaches; when using the continuous data, the unweighted estimate is substantially more precise.  




The R^2 estimates are shown in Table <ref>. As the heuristic analysis indicated, the ordinary pseudo-R^2 statistics are much larger than a design-based version, and so can be importantly biased relative to population or cohort statistics.  The bias exists both for the continuous-data models where the unweighted analysis is substantially more efficient, and the grouped-data models where there is little efficiency difference. 



 §.§ Informative two-phase sampling from a cohort



Wilms' Tumour is a rare, largely treatable childhood cancer of the kidney, and the National Wilms' Tumor Study Group has run a series of clinical trials aiming to reduce the long-term side-effects of treatment while maintaining the cure rate<cit.>. From a biostatistical viewpoint, one of the interesting aspects of Wilms' Tumour is that the histological classification (roughly, `cell abnormality') is difficult, and the study group central pathologist appears to be better at it than anyone else. That is, given the central lab histology, the local-hospital histology is not predictive of relapse, and the local-hospital histology can be treated as the central-lab histology plus random error.  There was obvious interest in studying Wilms' Tumour relapse with central-lab analysis of only a subset, rather than of all cases, and one of the data sets from the NWTG studies has become a standard example of two-phase sampling (eg, <cit.>)

Table <ref>, compares the results for the full cohort to two  sampling schemes:  case–control sampling with equal numbers of cases and controsl, and balanced two-phase sampling with the same number sampled from all four cells of a 2× 2 local-histology by relapse table.  The design-based estimator is always close to the full cohort estimator; however, the unweighted estimator under case–control sampling is biased upwards.  The bias is less dramatic than in the previous example because the control sampling fraction is larger: about 10% versus less than 1%. 






§ DISCUSSION


There is some controversy about the usefulness of pseudo-R^2 measures for generalised linear models, but they are quite widely used.  In this note I have shown how to construct design-based versions of the Nagelkerke and Cox–Snell pseudo-R^2, which should be of use to survey statisticians.  I have also shown that the standard versions of these statistics when used for logistic regression on case–control samples do not estimate the same model summary as they would under prospective or cross-sectional sampling — a fact that does not seem to be well known in biostatistics and epidemiology.  

abbrvnat



