SAS and R

In example 8.6 we showed how to change the reference category. This is the natural first thought analysts have when their primary comparisons aren't represented in the default output. But our interest might center on a number of comparisons which don't share a category. Or we might need to compare one group with the mean of the other groups. Today we'll explore tests and estimates of effects like these, which are sometimes called contrasts. In a later entry we'll explore more complex multivariate contrasts.

SAS

Access to this kind of comparison in SAS is provided in many model-fitting procedures using a test, estimate, or contrast statement. The differences among these can be subtle. In general, for simple comparisons, we recommend the estimate statement, where available. It is available for the important glm and genmod procedures, among others. In our example from linear regression, we changed the referent from heroin to alcohol by sorting the data and using the order=data option. This gains us pairwise comparisons between heroin and alcohol and between cocaine and alcohol. Here we show how to get the comparison between heroin and cocaine from the same model fit.

proc sort data=help_a; by descending substance; run;

proc glm data=help_a order=data;
 class substance;
 model i1 = age substance / solution;
 estimate "heroin vs. cocaine" substance -1 1 0 / e;
run; quit;

The syntax of the estimate statement is to optionally name the estimate (between the quotes) then to list the effects whose values we want to assess, followed by the desired values for each level of the effect. The e option requests that the contrast vector be printed-- this is a vital check that the contrast is working as desired. Here are the pieces of output generated by the estimate statement.

Coefficients for Estimate heroin vs. cocaine

 Row 1
Intercept 0

AGE 0

SUBSTANCE heroin -1
SUBSTANCE cocaine 1
SUBSTANCE alcohol 0

This is the result of the e option, and shows that we've correctly specified the difference between heroin and cocaine.

 Standard
Parameter Estimate Error t Value Pr> |t|
heroin vs. cocaine 3.00540049 2.15576257 1.39 0.1640

This is the comparison we want. It displays the difference, which we could confirm by examining the parameter estimates, as well as the standard error of the difference and the p-value, which can't be determined from the standard output.

To compare alcohol with the average of heroin and cocaine, we could use the following syntax (results omitted).

 estimate "cocaine + alcohol vs heroin" substance -2 1 1 /e divisor=2;

The divisor option allows us to use portions that are not easily represented as decimals. It's necessary because of the reserved use of the / in SAS. Any number of estimate statements can be issued in a single proc.

R

The fit.contrast() function in the gmodels package will generate contrasts for lm and glm objects. Multivariate contrasts can be found in the contrast() function in the Design package.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
lm3 = lm(i1 ~ substance + age, data=ds))
library(gmodels)
fit.contrast(lm3, "substance", c(0,-1, 1))

This generates the following output:

 Estimate Std. Error t value Pr(>|t|)
substance c=( 0 -1 1 ) -3.005400 2.155763 -1.394124 0.1639696

The simple syntax accepts model objects, the name of the effect, and a vector of contrasts to apply. The function does not replicate the contrast, which would be useful, but it is simple enough to check the parameter estimates from the model to ensure the desired result has been obtained.

The syntax for comparing heroin with the mean of alcohol and cocaine is straightforward.

> fit.contrast(lm3,"substance", c(-.5, -.5,1))
 Estimate Std. Error t value Pr(>|t|)
substance c=( -0.5 -0.5 1 ) -11.09965 1.904192 -5.829059 1.065763e-08 

The Hosmer and Lemeshow goodness of fit (GOF) test is a way to assess whether there is evidence for lack of fit in a logistic regression model. Simply put, the test compares the expected and observed number of events in bins defined by the predicted probability of the outcome. This can be calculated in R and SAS.

R

In R, we write a simple function to calculate the statistic and a p-value, based on vectors of observed and predicted probabilities. We use the cut() function (1.4.10) in concert with the quantile() function (2.1.5) to make the bins, then calculate the observed and expected counts, the chi-square statistic, and finally the associated p-value (Table 1.1). The function allows the user to define the number of bins but uses the common default of 10.

hosmerlem = function(y, yhat, g=10) {
 cutyhat = cut(yhat,
 breaks = quantile(yhat, probs=seq(0,
 1, 1/g)), include.lowest=TRUE)
 obs = xtabs(cbind(1 - y, y) ~ cutyhat)
 expect = xtabs(cbind(1 - yhat, yhat) ~ cutyhat)
 chisq = sum((obs - expect)^2/expect)
 P = 1 - pchisq(chisq, g - 2)
 return(list(chisq=chisq,p.value=P))
}

We'll run it with some of the HELP data (available at the book web site). Note that fitted(object) returns the predicted probabilities, if the object is the result of a call to glm() with family=binomial.

ds = read.csv("http://www.math.smith.edu/r/data/help.csv")
attach(ds)
logreg = glm(homeless ~ female + i1 + cesd + age + substance,
 family=binomial)
hosmerlem(y=homeless, yhat=fitted(logreg))

This returns the following output:

$chisq
[1] 8.495386
$p.value
[1] 0.3866328

The Design package, by Frank Harrell, includes the le Cessie and Houwelingen test (another goodness-of-fit test, Biometrics 1991 47:1267) and is also easy to run, though it requires using the package's function for logistic regression.

library(Design)
mod = lrm(homeless ~ female + i1 + cesd + age + substance, 
 x=TRUE, y=TRUE, data=ds)
resid(mod, 'gof')

Sum of squared errors Expected value|H0 SD
 104.1091804 103.9602955 0.1655883
 Z P
 0.8991269 0.3685851

The two tests are reassuringly similar.

SAS

In SAS, the Hosmer and Lemeshow goodness of fit test is generated with the lackfit option to the model statement in proc logistic (section 4.1.1). (We select out the results using the ODS system.)

ods select lackfitpartition lackfitchisq;
proc logistic data="c:\book\help.sas7bdat";
 class substance female;
 model homeless = female i1 cesd age substance / lackfit;
run;

This generates the following output:

 Partition for the Hosmer and Lemeshow Test

 HOMELESS = 1 HOMELESS = 0
 Group Total Observed Expected Observed Expected

 1 45 10 12.16 35 32.84
 2 45 12 14.60 33 30.40
 3 45 15 15.99 30 29.01
 4 45 17 17.20 28 27.80
 5 45 27 18.77 18 26.23
 6 45 20 20.28 25 24.72
 7 45 23 22.35 22 22.65
 8 45 25 25.04 20 19.96
 9 45 28 27.67 17 17.33
 10 48 32 34.95 16 13.05

Hosmer and Lemeshow Goodness-of-Fit Test

Chi-Square DF Pr> ChiSq

 8.4786 8 0.3882

The partition table shows the observed and expected count or events in each decile of the predicted probabilities.

The discrepancy between the SAS and R results is likely due to the odd binning SAS uses; the test is unstable in the presence of ties to the extent that some authorities suggest avoiding it. In general, with continuous predictors, the objections are not germane.