7 Statistical Analysis (SAS)

7.1 Setup

In this section, we’ll provide some examples in SAS of basic statistical analyses you might need to implement at some point.

The analyses you conduct as part of real research will almost certainly be more complicated than the code shown here. However, the hope is that you’ll get enough of a feel for the mechanics of SAS to begin putting the pieces together as demanded by a given project.

7.2 Exploratory Data Analysis

First, let’s read in our data and glance at it using a couple of different views.

SAS Documentation: PROC IMPORT

proc import datafile="data/nhefs.csv"
            out=nhefs
            dbms=csv
            replace;
run;

* information about the variables in NHEFS;
proc contents data=nhefs; run;

* print first five rows of dataset;
proc print data=nhefs(obs=5); run;

For now, let’s focus on observations (rows) that don’t have any missing data. To do this, let’s get a summary of how much missingness we have across all observations.

proc means data=nhefs nmiss STACKODSOUTPUT;
    ods output summary=nhefs_miss;
run;

proc print data=nhefs_miss; run;

Let’s drop every row with a missing value for any variable:

Complete case code adapted from an example on SAS blogs.

/* h/t https://blogs.sas.com/content/iml/2015/02/23/complete-cases.html */

data nhefs2;
    set nhefs;
    if cmiss(of _ALL_)=0;
run;

proc contents data=nhefs2; run;

Oops, we have only 44 rows left in the data set! Perhaps we should reconsider and pick several variables that will be of interest to us for our current purposes.

%let selectvars = age alcoholfreq asthma wt71 smokeintensity sbp allergies weakheart;

data nhefs3;
    set nhefs(keep=&selectvars);
    if cmiss(of _ALL_)=0;
run;

proc contents data=nhefs3; run;

That’s better. We’ve discarded 77 rows and focused on a subset of variables of interest.

Before we do some exploratory data analysis, let’s recall that category 5 for one of the variables we selected (alcoholfreq) indicated that alcohol frequency for that individual was unknown. This is missing data, and we should recode the alcoholfreq variable so that SAS recognizes these observations as such.

In addition, SAS currently believes that alcoholfreq is a numeric variable, though it is a categorical variable.

We can add a new variable to our nhefs3 data set that recodes alcoholfreq to meet our needs. We’ll name our new variable alcfreqcat

CHARACTER VARIABLES How did we get alcfreqcat to be a character variable when we based it on the numeric variable alcoholfreq? Notice in the data step that we initialized alcfreqcat using '.' instead of .. Wrapping the dot in quotes indicated to SAS that we wanted it to treat alfreqcat as a character variable.

data nhefs3;
    set nhefs3;
    if alcoholfreq = 5 then alcfreqcat = '.';
    else alcfreqcat = alcoholfreq;
run;

proc freq data=nhefs3;
    table alcfreqcat*alcoholfreq;
run;

proc contents data=nhefs3; run;

Although SAS will automatically throw out observations with missing data when running regressions and other procedures, it will do so based on which variables are included in model formulas. So that we conduct all the analyses below on the same subset of the NHEFS data, let’s drop any rows from our data set with missing values of alcfreqcat and drop the original alcoholfreq variable.

data nhefs4;
    set nhefs3(where=(alcfreqcat ne '.') drop=alcoholfreq);
run;

proc print data=nhefs4(obs=5); run;

Let’s summarize our dataset:

* save list of continuous variables and summarize;
%let contvars = age wt71 smokeintensity sbp;
proc means data=nhefs4;
    var &contvars;
run;

* save list of categorical variables and summarize;
%let catvars = asthma allergies alcfreqcat weakheart;
proc freq data=nhefs4;
    tables &catvars;
run;

7.3 Tabular Analysis

Soon you will be no stranger to two-by-two tables and contingency tables. Introductory statistics and epidemiologic methods courses will require you to analyze such tables, either to measure the strength of an exposure’s effect on some outcome or to establish a statistical association between two variables.

Going back to our NHEFS data, let’s say we’re interested in a possible association between quitting smoking and weight loss. The qsmk variable is a binary indicator of whether or not an individual quit smoking between 1971 and 1982, while the sbp variable measures an individual’s change in weight (kilograms) between 1971 and 1982.

Typically, we would not want to dichotomize a continuous variable like sbp, but for the sake of this example, let’s begin by creating a binary indicator of whether a subject gained weight during this timeframe.

data nhefs;
    set nhefs;
    if sbp > 140 then sbp_hi = 1;
    else if sbp <= 140 then sbp_hi = 0;
    if sbp = . then sbp_hi = .;
run;

proc freq data=nhefs;
    table sbp_hi;
run;

The summary statistics we generated above should comply with our expectations: those we marked as having a high systolic blood pressure have values of sbp above 140 mm/Hg, while those we marked as not having high systolic blood pressure all have values of sbp less than or equal to 140 mm/Hg.

Note, too, that we have 63 missing values in our new variable, just as we saw in the original outcome.

proc freq data=nhefs;
    table sbp_hi * qsmk /chisq;
run;

A natural choice here would be to conduct a two-sample test for the equality of proportions, where we compare the proportion of subjects in the quit smoking group with high systolic blood pressure versus the same proportion in the non-quit smoking group.

The output above shows us the group proportions, the \(\chi^2\) test statistic (elicited by specifying the table option chisq), and the corresponding p-value. We also get a 95% confidence interval for the difference in proportions.

7.4 Regression Models

In this section we will fit some regression models that you may encounter in classes or in your own research. We will leave the theoretical background to your statistics and methods courses but aim to provide enough information so that you understand the basic idea behind each approach and can refer back to this page in the future.

All upcoming examples use variables from the NHEFS dataset that we’ve referred back to several times by now.

ods graphics off;
proc reg data=nhefs;
    model sbp = qsmk /clb; * /clb adds confidence limits to output table;
quit;

data nhefs;
    set nhefs;
    sex_age = sex * age;
run;

ods graphics off;
proc reg data=nhefs;
    model sbp = qsmk smokeyrs smokeintensity diabetes sex age sex_age /clb;
quit;

proc genmod data=nhefs;
    model sbp = qsmk /link=identity dist=normal;
run;

proc genmod data=nhefs;
    model sbp = qsmk smokeyrs smokeintensity diabetes sex:age /link=identity dist=normal;
run;

proc logistic data=nhefs;
    model sbp_hi(event='1') = qsmk;
    estimate 'qsmk (1 vs. 0)' qsmk 1 /exp;
run;

proc genmod data=nhefs;
    model sbp_hi(ref='0') = qsmk /link=logit dist=bin;
    estimate 'qsmk (1 vs. 0)' qsmk 1 /exp;
run;

proc genmod data=nhefs;
    model sbp_hi(ref='0') = qsmk /link=log dist=bin;
    estimate 'qsmk (1 vs. 0)' qsmk 1 /exp;
run;

proc genmod data=nhefs;
    class seqn;
    model sbp_hi(ref='0') = qsmk /link=log dist=poisson;
    repeated subject=seqn;
    estimate 'qsmk (1 vs. 0)' qsmk 1 /exp;
run;