Regression with SAS
Chapter 4 - Beyond OLS

Chapter Outline
    4.1 Robust Regression Methods
        4.1.1 Regression with Robust Standard Errors
        4.1.2 Using the Proc Genmod for Clustered Data
        4.1.3 Robust Regression
        4.1.4 Quantile Regression
    4.2 Constrained Linear Regression
    4.3 Regression with Censored or Truncated Data
        4.3.1 Regression with Censored Data
        4.3.2 Regression with Truncated Data
    4.4 Regression with Measurement Error
    4.5 Multiple Equation Regression Models
        4.5.1 Seemingly Unrelated Regression
        4.5.2 Multivariate Regression
    4.6 Summary     

In this chapter we will go into various commands that go beyond OLS. This chapter is a bit different from the others in that it covers a number of different concepts, some of which may be new to you. These extensions, beyond OLS, have much of the look and feel of OLS but will provide you with additional tools to work with linear models.

The topics will include robust regression methods, constrained linear regression, regression with censored and truncated data, regression with measurement error, and multiple equation models.

4.1 Robust Regression Methods

It seems to be a rare dataset that meets all of the assumptions underlying multiple regression. We know that failure to meet assumptions can lead to biased estimates of coefficients and especially biased estimates of the standard errors. This fact explains a lot of the activity in the development of robust regression methods.

The idea behind robust regression methods is to make adjustments in the estimates that take into account some of the flaws in the data itself. We are going to look at three robust methods: regression with robust standard errors, regression with clustered data, robust regression, and quantile regression.

Before we look at these approaches, let's look at a standard OLS regression using the elementary school academic performance index (elemapi2.dta) dataset. We will look at a model that predicts the api 2000 scores using the average class size in K through 3 (acs_k3), average class size 4 through 6 (acs_46), the percent of fully credentialed teachers (full), and the size of the school (enroll). First let's look at the descriptive statistics for these variables.  Note the missing values for acs_k3 and acs_k6.

 proc means data = "c:sasregelemapi2" mean std max min;   var api00 acs_k3 acs_46 full enroll; run;

The MEANS ProcedureVariable            Mean         Std Dev         Minimum         Maximum------------------------------------------------------------------------api00        647.6225000     142.2489610     369.0000000     940.0000000acs_k3        19.1608040       1.3686933      14.0000000      25.0000000acs_46        29.6851385       3.8407840      20.0000000      50.0000000full          84.5500000      14.9497907      37.0000000     100.0000000enroll       483.4650000     226.4483847     130.0000000         1570.00------------------------------------------------------------------------

Below we see the regression predicting api00 from acs_k3 acs_46 full and enroll. We see that all of the variables are significant except for acs_k3.

 proc reg data = "c:sasregelemapi2";   model api00 = acs_k3 acs_46 full enroll ;  run;

The REG ProcedureModel: MODEL1Dependent Variable: api00                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     4        3071909         767977      61.01    <.0001Error                   390        4909501          12588Corrected Total         394        7981410Root MSE            112.19832    R-Square     0.3849Dependent Mean      648.65063    Adj R-Sq     0.3786Coeff Var            17.29719                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1       -5.20041       84.95492      -0.06      0.9512acs_k3        1        6.95438        4.37110       1.59      0.1124acs_46        1        5.96601        1.53105       3.90      0.0001full          1        4.66822        0.41425      11.27      <.0001enroll        1       -0.10599        0.02695      -3.93      <.0001

Since the regression procedure is interactive and we haven't issued the quit command, we can test both of the class size variables, and we find the overall test of these two variables is significant.

test acs_k3 = acs_46 = 0; run;

       Test 1 Results for Dependent Variable api00                                MeanSource             DF         Square    F Value    Pr > FNumerator           2         139437      11.08    <.0001Denominator       390          12588

Here is the residual versus fitted plot for this regression. Notice that the pattern of the residuals is not exactly as we would hope.  The spread of the residuals is somewhat wider toward the middle right of the graph than at the left, where the variability of the residuals is somewhat smaller, suggesting some heteroscedasticity.

plot r.*p.; run;  

Here is the index plot of Cook's D for this regression. We see 4 points that are somewhat high in both their leverage and their residuals. 

  plot cookd.*obs.; run;  

None of these results are dramatic problems, but the plot of residual vs. predicted value suggests that there might be some outliers and some possible heteroscedasticity and the index plot of Cook's D shows some points in the upper right quadrant that could be influential. We might wish to use something other than OLS regression to estimate this model. In the next several sections we will look at some robust regression methods.

4.1.1 Regression with Robust Standard Errors

The SAS proc reg  includes an option called acov in the model statement for estimating the asymptotic covariance matrix of the estimates under the hypothesis of heteroscedasticity. The standard error obtained from the asymptotic covariance matrix is considered to be more  robust and can deal with a collection of minor concerns about failure to meet assumptions, such as minor problems about normality, heteroscedasticity, or some observations that exhibit large residuals, leverage or influence. For such minor problems, the standard error based on acov  may effectively deal with these concerns.

With the acov option, the point estimates of the coefficients are exactly the same as in ordinary OLS, but we will calculate the standard errors based on the asymptotic covariance matrix. Here is the same regression as above using the acov option. We also use SAS ODS (Output Delivery System)  to output the parameter estimates along with the asymptotic covariance matrix. We calculated the robust standard error in a data step and merged them with the parameter estimate using proc sql and created the t-values and corresponding probabilities.  Note the changes in the standard errors and t-tests (but no change in the coefficients). In this particular example, using robust standard errors did not change any of the conclusions from the original OLS regression. We should also mention that the robust standard error has been adjusted for the sample size correction. 

proc reg data = "c:sasregelemapi2";       model api00 = acs_k3 acs_46 full enroll /acov;       ods output  ACovEst = estcov;       ods output ParameterEstimates=pest; run; quit; data temp_dm;   set estcov;   drop model dependent;   array a(5) intercept acs_k3 acs_46 full enroll;   array b(5) std1-std5;   b(_n_) = sqrt((395/390)*a(_n_));   std = max(of std1-std5);   keep variable std; run; proc sql;   select pest.variable, estimate, stderr, tvalue, probt, std as robust_stderr,           estimate/robust_stderr as tvalue_rb,           (1 - probt(abs(estimate/robust_stderr), 394))*2 as probt_rb    from pest, temp_dm   where pest.variable=temp_dm.variable; quit;

                                                      robust_Variable   Estimate     StdErr  tValue   Probt     stderr  tvalue_rb  probt_rb------------------------------------------------------------------------------Intercept  -5.20041   84.95492   -0.06  0.9512   86.66308   -0.06001   0.95218acs_k3      6.95438    4.37110    1.59  0.1124   4.620599   1.505082  0.133104acs_46      5.96601    1.53105    3.90  0.0001   1.573214   3.792246  0.000173full        4.66822    0.41425   11.27  <.0001   0.414681   11.25737         0enroll     -0.10599    0.02695   -3.93  <.0001   0.028015   -3.78331  0.000179

4.1.2 Using the Proc Genmod for Clustered Data

As described in Chapter 2, OLS regression assumes that the residuals are independent. The elemapi2 dataset contains data on 400 schools that come from 37 school districts. It is very possible that the scores within each school district may not be independent, and this could lead to residuals that are not independent within districts. SAS proc genmod is used to model correlated data.  We can use the class statement and the repeated statement  to indicate that the observations are clustered into districts (based on dnum) and that the observations may be correlated within districts, but would be independent between districts. 

proc genmod data="c:sasregelemapi2";     class dnum;     model api00 = acs_k3 acs_46 full enroll ;     repeated subject=dnum / type=ind ; run; quit;

The GENMOD Procedure

             Analysis Of GEE Parameter Estimates              Empirical Standard Error Estimates

                   Standard   95% ConfidenceParameter Estimate    Error       Limits            Z Pr > |Z|

Intercept  -5.2004 119.5172 -239.450 229.0490   -0.04   0.9653acs_k3      6.9544   6.7726  -6.3196  20.2284    1.03   0.3045acs_46      5.9660   2.4839   1.0976  10.8344    2.40   0.0163full        4.6682   0.6904   3.3151   6.0213    6.76   <.0001enroll     -0.1060   0.0421  -0.1886  -0.0234   -2.51   0.0119

As with the regression with robust error, the estimate of the coefficients are the same as the OLS estimates, but the standard errors take into account that the observations within districts are non-independent.  Even though the standard errors are larger in this analysis, the three variables that were significant in the OLS analysis are significant in this analysis as well. 

We notice that the standard error estimates given here are different from what Stata's result using regress with the cluster option. This is because that Stata further does a finite-sample adjustment. We can do some SAS programming here for the adjustment. The adjusted variance is a constant times the variance obtained from the empirical standard error estimates. This particular constant is
(N-1)/(N-k)*M/(M-1).

data em;   set 'c:sasregelemapi2'; run;

proc genmod data=em;     class dnum;     model api00 = acs_k3 acs_46 full enroll ;     repeated subject=dnum / type = ind covb ; 	ods output geercov = gcov;     ods output GEEEmpPEst = parms; run; quit;

proc sql;   select count(dnum),count(distinct dnum)  into :n, :m   from em; quit; proc sql;   select count(prm1) into :k   from gcov; quit; data gcov_ad;   set gcov;   array all(*) _numeric_;   do i = 1 to dim(all);   all(i) = all(i)*((&n-1)/(&n-&k))*(&m/(&m-1));   if i = _n_ then std_ad = sqrt(all(i));   end; drop i; keep std_ad; run; data all;   merge parms gcov_ad; run; proc print data = all noobs; run;

Parm         Estimate      Stderr     LowerCL     UpperCL          Z     ProbZ     std_ad

Intercept     -5.2004    119.5172    -239.450    229.0490      -0.04    0.9653    121.778acs_k3         6.9544      6.7726     -6.3196     20.2284       1.03    0.3045      6.901acs_46         5.9660      2.4839      1.0976     10.8344       2.40    0.0163      2.531full           4.6682      0.6904      3.3151      6.0213       6.76    <.0001      0.703enroll        -0.1060      0.0421     -0.1886     -0.0234      -2.51    0.0119      0.043

4.1.3 Robust Regression

In SAS, we can not simply execute some proc to perform a robust regression using iteratively reweighted least squares. In order to perform a robust regression,  we have to write our own macro. To this end, ATS has written a macro called robust_hb.sas. This macro first uses Hubert weight and later switches to biweight. This is why the macro is called robust_hb where h and b stands for Hubert and biweight respectively.  Robust regression assigns a weight to each observation with higher weights given to better behaved observations. In fact, extremely deviant cases, those with Cook's D greater than 1, can have their weights set to missing so that they are not included in the analysis at all.

The macro robust_hb.sas uses another macro called mad.sas to generate MAD (median absolute deviation) during the iteration process. We will include both macros to perform the robust regression analysis as shown below. Note that in this analysis both the coefficients and the standard errors differ from the original OLS regression. 

%include 'c:sasregmad.sas'; %include 'c:sasregrobust_hb.sas'; %robust_hb("c:sasregelemapi2", api00, acs_k3 acs_46 full enroll,             .01, 0.00005, 10);

The REG ProcedureModel: MODEL1Dependent Variable: api00Weight: _w2_                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     4        3053248         763312      73.19    <.0001Error                   390        4067269          10429Corrected Total         394        7120516Root MSE            102.12196    R-Square     0.4288Dependent Mean      647.38090    Adj R-Sq     0.4229Coeff Var            15.77463                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1       -6.79858       80.45406      -0.08      0.9327acs_k3        1        6.11031        4.15015       1.47      0.1417acs_46        1        6.25588        1.45280       4.31      <.0001full          1        4.79575        0.39212      12.23      <.0001enroll        1       -0.10924        0.02558      -4.27      <.0001

If you compare the robust regression results (directly above) with the OLS results previously presented, you can see that the coefficients and standard errors are quite similar, and the t values and p values are also quite similar. Despite the minor problems that we found in the data when we performed the OLS analysis, the robust regression analysis yielded quite similar results suggesting that indeed these were minor problems.   Had the results been substantially different, we would have wanted to further investigate the reasons why the OLS and robust regression results were different, and among the two results the robust regression results would probably be the more trustworthy.

Let's look at the predicted (fitted) values (p), the residuals (r), and the leverage (hat) values (h).  The macro robust_hb generates a final data set with predicted values, raw residuals and leverage values together with the original data called _tempout_.

Now, let's check on the various predicted values and the weighting. First, we will sort by _w2_, the weight generated at last iteration. Then we will look at the first 15 observations. Notice that the smallest weights are near one-half but quickly get into the .6 range. The first five values are missing due to the missing values of predictors. 

 proc sort data = _tempout_;   by descending _w2_; run; proc print data = _tempout_ (obs=10);   var snum api00 p r h _w2_; run;

Obs    snum    api00       p           r           h         _w2_  1     116     513        .           .        .           .  2    3072     763        .           .        .           .  3    3055     590        .           .        .           .  4    4488     521        .           .        .           .  5    4534     445        .           .        .           .  6     637     447     733.426    -286.144    0.003657    0.55684  7    5387     892     611.709     280.464    0.002278    0.57183  8    2267     897     622.062     275.512    0.009887    0.58481  9      65     903     631.930     271.731    0.010198    0.59466 10    3759     585     842.664    -257.473    0.040009    0.63128

Now, let's look at the last 10 observations. The weights for observations with snum 1678, 4486 and 1885 are all very close to one, since the residuals are fairly small. Note that the observations above that have the lowest weights are also those with the largest residuals (residuals over 200) and the observations below with the highest weights have very low residuals (all less than 3).

proc sort data = _tempout_;   by _w2_; run; proc print data = _tempout_ (obs=10);   var snum api00 p r h _w2_; run;

Obs    snum    api00       p           r           h         _w2_  1    1678     497     497.065     0.18424    0.023137    1.00000  2    4486     706     706.221     0.20151    0.013740    1.00000  3    1885     605     606.066    -0.41421    0.013745    1.00000  4    3535     705     703.929     1.16143    0.004634    0.99999  5    3024     727     729.278    -2.01559    0.010113    0.99997  6    3700     717     719.803    -2.43802    0.007317    0.99996  7    1596     536     533.918     2.77985    0.012143    0.99995  8     181     724     728.068    -3.53209    0.013665    0.99992  9    1949     696     691.898     3.96736    0.020426    0.99990 10    5917     735     731.141     4.03336    0.005831    0.99990

After using macro robust_hb.sas, we can use the dataset _tempout_ to create some graphs for regression diagnostic purposes. For example,  we can create a graph of residuals versus fitted (predicted) with a line at zero. This plot looks much like the OLS plot, except that in the OLS all of the observations would be weighted equally, but as we saw above the observations with the greatest residuals are weighted less and hence have less influence on the results.

axis1 order = (-300 to 300 by 100) label=(a=90) minor=none; axis2 order = (300 to 900 by 300) minor=none; symbol v=star h=0.8 c=blue; proc gplot data = _tempout_;   plot r*p = 1 /haxis=axis2 vaxis=axis1 vref=0;   label r ="Residual";   label p ="Fitted Value"; run; quit;

4.1.4 Quantile Regression

Quantile regression, in general, and median regression, in particular, might be considered as an alternative to robust regression. SAS does quantile regression using a little bit of proc iml.  Inside proc iml, a procedure called LAV is called and it does a median regression in which the coefficients will be estimated by minimizing the absolute deviations from the median. Of course, as an estimate of central tendency, the median is a resistant measure that is not as greatly affected by outliers as is the mean. It is not clear that median regression is a resistant estimation procedure, in fact, there is some evidence that it can be affected by high leverage values.

Here is what the quantile regression looks like using SAS proc iml. The first data step is to make sure that the data set that proc iml takes as input does not have any missing values. Inside proc iml we first generate necessary matrices for regression computation and then call the procedure LAV. After calling LAV we can calculate the predicted values and residuals. At last, we create a data set called _temp_ containing the dependent variables and all the predictors plus the predicted values and residuals. 

data elemapi2;    set "c:sasregelemapi2";    cons = 1;    if api00 ~=. & acs_k3 ~= . & acs_46 ~=. & full ~=. & enroll ~=.; run;  proc iml ; /*Least absolute values*/   use elemapi2;   read all;   a = cons || acs_k3 || acs_46 || full || enroll;   b=api00;   opt= { . 3  0 1 };    call lav(rc,xr,a,b,,opt); /*print out the estimates*/      opt= { . 0  . 1 };    call lav(rc,xr,a,b,,opt); /*no print out but to create the xr*/     pred = a*t(xr);   resid = b - pred;   create _temp_ var { api00 cons acs_k3 acs_46 full enroll pred resid};   append; quit;

                                 LAV (L1) Estimation                               Start with LS Solution                      Start Iter: gamma=284.75134813 ActEqn=395 Iter  N Huber  Act Eqn     Rank    Gamma         L1(x)      F(Gamma)   31       36        5        5   0.0000    36268.1049    36240.6335Algorithm convergedObjective Function L1(x)= 36268.104941Number Search Directions=       68Number Refactorizations =        2Total Execution Time=       0.0200Necessary and sufficient optimality condition satisfied.                                    L1 Solution with ASEEst  17.1505029  1.2690656888   7.2240793844   5.3238408715 -0.124573396ASE  123.531545  6.3559394265   2.2262729207   0.602359504   0.0391932684

The coefficient and standard error for acs_k3 are considerably different as compared to OLS (the coefficients are 1.2 vs 6.9 and the standard errors are 6.4 vs 4.3). The coefficients and standard errors for the other variables are also different, but not as dramatically different. Nevertheless, the  quantile regression results indicate that, like the OLS results, all of the variables except acs_k3 are significant. Using the data set _temp_ we created above we obtain a plot of residuals vs. predicted values shown below.

symbol v=star h=0.8 c=blue; axis1 order = (-300 to 300 by 100) label=(a=90) minor=none; axis2 order = (300 to 900 by 300) minor=none; proc gplot data = _temp_;   plot resid*pred = 1 / vaxis=axis1 haxis=axis2 vref=0;   label resid="Residual";   label pred = "Fitted Value"; run; quit;

4.2 Constrained Linear Regression

Let's begin this section by looking at a regression model using the hsb2 dataset.  The hsb2 file is a sample of 200 cases from the Highschool and Beyond Study (Rock, Hilton, Pollack, Ekstrom & Goertz, 1985). It includes the following variables: id female race ses schtyp program read write math science socst. The variables read write math science socst are the results of standardized tests on reading, writing, math, science and social studies (respectively), and the variable female is coded 1 if female, 0 if male.

Let's start by doing an OLS regression where we predict socst score from read, write, math, science and female (gender)

proc reg data="c:sasreghsb2";   model  socst = read write math science female ; run;

The REG ProcedureModel: MODEL1Dependent Variable: socst                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     5          10949     2189.85150      35.44    <.0001Error                   194          11987       61.78834Corrected Total         199          22936Root MSE              7.86056    R-Square     0.4774Dependent Mean       52.40500    Adj R-Sq     0.4639Coeff Var            14.99963                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1        7.33934        3.65024       2.01      0.0457read          1        0.37840        0.08063       4.69      <.0001write         1        0.38587        0.08893       4.34      <.0001math          1        0.13033        0.08938       1.46      0.1464science       1       -0.03339        0.08187      -0.41      0.6838female        1       -0.35326        1.24537      -0.28      0.7770

Notice that the coefficients for read and write are very similar, which makes sense since they are both measures of language ability.  Also, the coefficients for math and science are similar (in that they are both not significantly different from 0). Suppose that we have a theory that suggests that read and write and math should have equal coefficients. We can test the equality of the coefficients using the test command.

   test read = write;  run;

        Test 1 Results for Dependent Variable socst                                MeanSource             DF         Square    F Value    Pr > FNumerator           1        0.19057       0.00    0.9558Denominator       194       61.78834

The test result indicates that there is no significant difference in the coefficients for the reading and writing scores. Since it appears that the coefficients for math and science are also equal, let's test the equality of those as well.

   test math = science;   run;

       Test 2 Results for Dependent Variable socst                                MeanSource             DF         Square    F Value    Pr > FNumerator           1       89.63950       1.45    0.2299Denominator       194       61.78834

Let's now perform both of these tests together, simultaneously testing that the coefficient for read equals write and math equals science using mtest statement.  

  mtest math - science, read - write;   run;

Multivariate Test 1                Multivariate Statistics and Exact F Statistics                              S=1    M=0    N=96Statistic                        Value    F Value    Num DF    Den DF    Pr > FWilks' Lambda               0.99257173       0.73         2       194    0.4852Pillai's Trace              0.00742827       0.73         2       194    0.4852Hotelling-Lawley Trace      0.00748386       0.73         2       194    0.4852Roy's Greatest Root         0.00748386       0.73         2       194    0.4852

Note this second test has 2 df, since it is testing both of the hypotheses listed, and this test is not significant, suggesting these pairs of coefficients are not significantly different from each other.  We can estimate regression models where we constrain coefficients to be equal to each other.  For example, let's begin on a limited scale and constrain read to equal write. Proc reg uses restrict statement to  accomplish this.

proc reg data="c:sasreghsb2";   model  socst = read write math science female ;   restrict read=write; run;

The REG ProcedureModel: MODEL1Dependent Variable: socstNOTE: Restrictions have been applied to parameter estimates.                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     4          10949     2737.26674      44.53    <.0001Error                   195          11987       61.47245Corrected Total         199          22936Root MSE              7.84044    R-Square     0.4774Dependent Mean       52.40500    Adj R-Sq     0.4667Coeff Var            14.96124                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1        7.35415        3.63118       2.03     0.0442read          1        0.38185        0.05139       7.43     <.0001write         1        0.38185        0.05139       7.43     <.0001math          1        0.13030        0.08915       1.46     0.1454science       1       -0.03328        0.08164      -0.41     0.6840female        1       -0.32962        1.16736      -0.28     0.7780RESTRICT     -1      -25.51302      458.21611      -0.06     0.9558** Probability computed using beta distribution.

Notice that the coefficients for read and write are identical, along with their standard errors, t-test, etc. Also note that the degrees of freedom for the F test is four, not five, as in the OLS model. This is because only one coefficient is estimated for read and write, estimated like a single variable equal to the sum of their values.  Notice also that the Root MSE is slightly higher for the constrained model, but only slightly higher.  This is because we have forced the model to estimate the coefficients for read and write that are not as good at minimizing the Sum of Squares Error (the coefficients that would minimize the SSE would be the coefficients from the unconstrained model).

Next, we will define a second constraint, setting math equal to science together with the first constraint we set before. 

proc reg data="c:sasreghsb2";   model  socst = read write math science female ;   restrict read = write, math = science; run;

The REG ProcedureModel: MODEL1Dependent Variable: socstNOTE: Restrictions have been applied to parameter estimates.                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     3          10860     3619.84965      58.75    <.0001Error                   196          12077       61.61554Corrected Total         199          22936Root MSE              7.84956    R-Square     0.4735Dependent Mean       52.40500    Adj R-Sq     0.4654Coeff Var            14.97864                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1        7.50566        3.63323       2.07     0.0402read          1        0.38604        0.05133       7.52     <.0001write         1        0.38604        0.05133       7.52     <.0001math          1        0.04281        0.05192       0.82     0.4107science       1        0.04281        0.05192       0.82     0.4107female        1       -0.20087        1.16383      -0.17     0.8631RESTRICT     -1      -15.36285      458.82638      -0.03     0.9734*RESTRICT     -1      547.24384      454.01563       1.21     0.2290** Probability computed using beta distribution.

Now the coefficients for read =  write and math = science and the degrees of freedom for the model has dropped to three.  Again, the Root MSE is slightly larger than in the prior model, but we should emphasize only very slightly larger.  If indeed the population coefficients for read =  write and math = science, then these combined (constrained) estimates may be more stable and generalize better to other samples.  So although these estimates may lead to slightly higher standard error of prediction in this sample, they may generalize better to the population from which they came.

4.3 Regression with Censored or Truncated Data

Analyzing data that contain censored values or are truncated is common in many research disciplines. According to Hosmer and Lemeshow (1999), a censored value is one whose value is incomplete due to random factors for each subject. A truncated observation, on the other hand, is one which is incomplete due to a selection process in the design of the study.

We will begin by looking at analyzing data with censored values.

4.3.1 Regression with Censored Data

In this example we have a variable called acadindx which is a weighted combination of standardized test scores and academic grades. The maximum possible score on acadindx is 200 but it is clear that the 16 students who scored 200 are not exactly equal in their academic abilities. In other words, there is variability in academic ability that is not being accounted for when students score 200 on acadindx. The variable acadindx is said to be censored, in particular, it is right censored.

Let's look at the example. We will begin by looking at a description of the data, some descriptive statistics, and correlations among the variables.

proc means data = "c:sasregacadindx"; run;

The MEANS ProcedureVariable      N            Mean         Std Dev         Minimum         Maximum-------------------------------------------------------------------------------id          200     100.5000000      57.8791845       1.0000000     200.0000000female      200       0.5450000       0.4992205               0       1.0000000reading     200      52.2300000      10.2529368      28.0000000      76.0000000writing     200      52.7750000       9.4785860      31.0000000      67.0000000acadindx    200     172.1850000      16.8173987     138.0000000     200.0000000-------------------------------------------------------------------------------

proc means data = "c:sasregacadindx" N;   where acadindx = 200; run;

The MEANS ProcedureVariable      N---------------id           16female       16reading      16writing      16acadindx     16---------------

proc corr data = "c:sasregacadindx" nosimple noprob;   var acadindx female reading writing; run;

The CORR Procedure   4  Variables:    acadindx female   reading  writing           Pearson Correlation Coefficients, N = 200              acadindx        female       reading       writingacadindx       1.00000      -0.08210       0.71309       0.66256female        -0.08210       1.00000      -0.05308       0.25649reading        0.71309      -0.05308       1.00000       0.59678writing        0.66256       0.25649       0.59678       1.00000

Now, let's run a standard OLS regression on the data and generate predicted scores in p1.

proc reg data = "c:sasregacadindx";    model acadindx =female reading writing;    output out = reg1 p = p1; run; quit;

The REG ProcedureModel: MODEL1Dependent Variable: acadindx                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     3          34994          11665     107.40    <.0001Error                   196          21288      108.61160Corrected Total         199          56282Root MSE             10.42169    R-Square     0.6218Dependent Mean      172.18500    Adj R-Sq     0.6160Coeff Var             6.05261                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1       96.11841        4.48956      21.41      <.0001female        1       -5.83250        1.58821      -3.67      0.0003reading       1        0.71842        0.09315       7.71      <.0001writing       1        0.79057        0.10410       7.59      <.0001

The proc lifereg is one of the procedures in SAS that can be used for regression with censored data. The syntax of the command is similar to proc reg with the addition of the variable indicating if an observation is censored. Therefore, we have to create a data set with the information on censoring. 

data tobit_model;   set "c:sasregacadindx";   censor = ( acadindx >= 200 ); run;   proc lifereg data = tobit_model;    model acadindx*censor(1) = female reading writing /d=normal;    output out = reg2 p = p2; run;

The LIFEREG Procedure             Model InformationData Set                    WORK.TOBIT_MODELDependent Variable                  acadindxCensoring Variable                    censorCensoring Value(s)                         1Number of Observations                   200Noncensored Values                       184Right Censored Values                     16Left Censored Values                       0Interval Censored Values                   0Name of Distribution                  NormalLog Likelihood                  -718.0636168Algorithm converged.       Type III Analysis of Effects                         WaldEffect       DF    Chi-Square    Pr > ChiSqfemale        1       14.0654        0.0002reading       1       60.8529        <.0001writing       1       54.1655        <.0001                  Analysis of Parameter Estimates                      Standard   95% Confidence     Chi-Parameter DF Estimate    Error       Limits       Square Pr > ChiSqIntercept  1  92.7378   4.8034  83.3233 102.1524  372.74     <.0001female     1  -6.3473   1.6924  -9.6644  -3.0302   14.07     0.0002reading    1   0.7777   0.0997   0.5823   0.9731   60.85     <.0001writing    1   0.8111   0.1102   0.5951   1.0271   54.17     <.0001Scale      1  10.9897   0.5817   9.9067  12.1912

Let's merge the two data sets we created together to compare the predicted score p1 and p2.   

data compare;   merge reg1 reg2;   by id; run; proc means data = compare;   var acadindx p1 p2; run;

The MEANS ProcedureVariable      N            Mean         Std Dev         Minimum         Maximum-------------------------------------------------------------------------------acadindx    200     172.1850000      16.8173987     138.0000000     200.0000000p1          200     172.1850000      13.2608696     142.3820725     201.5311070p2          200     172.7040275      14.0029221     141.2210940     203.8540576-------------------------------------------------------------------------------

It shows that the censored regression model predicted values have a larger standard deviation and a greater range of values. When we look at a listing of p1 and p2 for all students who scored the maximum of 200 on acadindx, we see that in every case the censored regression model predicted value is greater than the OLS predicted value. These predictions represent an estimate of what the variability would be if the values of acadindx could exceed 200.

proc print data = compare;   var acadindx p1 p2;   where acadindx = 200; run;

Obs    acadindx       p1         p2 32       200      179.175    179.620 57       200      192.681    194.329 68       200      201.531    203.854 80       200      191.831    193.577 82       200      188.154    189.563 88       200      186.573    187.940 95       200      195.997    198.176100       200      186.933    188.108132       200      197.578    199.798136       200      189.459    191.144143       200      191.185    192.833157       200      191.614    193.477161       200      180.251    181.008169       200      182.275    183.367174       200      191.614    193.477200       200      187.662    189.421

Proc lifereg  handles right censored, left censored and interval censored data.  There are two other commands in SAS that perform censored regression analysis such as proc qlim. 

4.3.2 Regression with Truncated Data

Truncated data occurs when some observations are not included in the analysis because of the value of the variable. We will illustrate analysis with truncation using the dataset, acadindx, that was used in the previous section. 

Let's imagine that in order to get into a special honors program, students need to score at least 160 on acadindx. So we will drop all observations in which the value of acadindx is less than or equal 160. Now, let's estimate the same model that we used in the section on censored data, only this time we will pretend that a 200 for acadindx is not censored.

proc reg data = "c:sasregacadindx";   model acadindx = female reading writing;   where acadindx >160; run; quit;

The REG ProcedureModel: MODEL1Dependent Variable: acadindx                             Analysis of Variance                                    Sum of           MeanSource                   DF        Squares         Square    F Value    Pr > FModel                     3     8074.79638     2691.59879      33.01    <.0001Error                   140          11416       81.54545Corrected Total         143          19491Root MSE              9.03025    R-Square     0.4143Dependent Mean      180.42361    Adj R-Sq     0.4017Coeff Var             5.00503                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1      125.63547        5.89156      21.32      <.0001female        1       -5.23849        1.61563      -3.24      0.0015reading       1        0.44111        0.09635       4.58      <.0001writing       1        0.58733        0.11508       5.10      <.0001

It is clear that the estimates of the coefficients are distorted due to the fact that 53 observations are no longer in the dataset. This amounts to restriction of range on both the response variable and the predictor variables. For example, the coefficient for writing dropped from .79 to .58. What this means is that if our goal is to find the relation between adadindx and the predictor variables in the populations, then the truncation of acadindx in our sample is going to lead to biased estimates. A better approach to analyzing these data is to use truncated regression. In SAS this can be accomplished using proc qlim. Proc qlim is an experimental procedure first available in SAS version 8.1. Proc qlim (Qualitative and LImited dependent variable model) analyzes univariate (and multivariate) limited dependent variable models where dependent variables takes discrete values or dependent variables are observed only in a limited range of values. 

data trunc_model;   set "c:sasregacadindx";   y = .;   if acadindx > 160 & acadindx ~=. then y = acadindx; run; proc qlim data=trunc_model ;   model y  = female reading writing;   endogenous y~ truncated (lb=160); run;

The QLIM Procedure                            Summary Statistics of Continuous Responses                                                                                    N Obs    N Obs                            Standard                           Lower       Upper    Lower    UpperVariable        Mean           Error          Type             Bound       Bound    Bound    Bound   y        180.4236       11.674837       Truncated             160             Model Fit SummaryNumber of Endogenous Variables             1Endogenous Variable                        yNumber of Observations                   144Missing Values                            56Log Likelihood                    -510.00768Maximum Absolute Gradient         3.87471E-6Number of Iterations                      14AIC                                     1030Schwarz Criterion                       1045Algorithm converged.                      Parameter Estimates                                 Standard                 ApproxParameter        Estimate           Error    t Value    Pr > |t|Intercept      110.289206        8.673847      12.72      <.0001female          -6.099602        1.925245      -3.17      0.0015reading          0.518179        0.116829       4.44      <.0001writing          0.766164        0.152620       5.02      <.0001_Sigma           9.803571        0.721646      13.59      <.0001

The coefficients from the proc qlim are closer to the OLS results, for example the coefficient for writing is .77 which is closer to the OLS results of .79.  However, the results are still somewhat different on the other variables, for example the coefficient for reading is .52 in the proc qlim as compared to .72 in the original OLS with the unrestricted data, and better than the OLS estimate of .47 with the restricted data.  While proc qlim may improve the estimates on a restricted data file as compared to OLS, it is certainly no substitute for analyzing the complete unrestricted data file.

4.4 Regression with Measurement Error

As you will most likely recall, one of the assumptions of regression is that the predictor variables are measured without error. The problem is that measurement error in predictor variables leads to under estimation of the regression coefficients. 

This section is under development.

4.5 Multiple Equation Regression Models

If a dataset has enough variables we may want to estimate more than one regression model. For example, we may want to predict y1 from x1 and also predict y2 from x2. Even though there are no variables in common these two models are not independent of one another because the data come from the same subjects. This is an example of one type multiple equation regression known as seemly unrelated regression.. We can estimate the coefficients and obtain standard errors taking into account the correlated errors in the two models. An important feature of multiple equation modes is that we can test predictors across equations.

Another example of multiple equation regression is if we wished to predict y1, y2 and y3 from x1 and x2. This is a three equation system, known as multivariate regression, with the same predictor variables for each model. Again, we have the capability of testing coefficients across the different equations.

Multiple equation models are a powerful extension to our data analysis tool kit.

4.5.1 Seemingly Unrelated Regression

Let's continue using the hsb2 data file to illustrate the use of seemingly unrelated regression. This time let's look at two regression models.

    science = math female     write   = read female    

It is the case that the errors (residuals) from these two models would be correlated. This would be true even if the predictor female were not found in both models. The errors would be correlated because all of the values of the variables are collected on the same set of observations. This is a situation tailor made for seemingly unrelated regression using the proc syslin with option sur. With the proc syslin we can estimate both models simultaneously while accounting for the correlated errors at the same time, leading to efficient estimates of the coefficients and standard errors. The syntax is as follows.

proc syslin data="c:sasreghsb2" sur ;   science:  model science = math female ;   write:    model write = read female ; run;

The first part of the output consists of the OLS estimate for each model. Here is the OLS estimate for the first model. 

               The SYSLIN ProcedureOrdinary Least Squares EstimationModel                  SCIENCEDependent Variable     science                         Analysis of Variance                               Sum of        MeanSource                 DF     Squares      Square    F Value    Pr > FModel                   2    7993.550    3996.775      68.38    <.0001Error                 197    11513.95    58.44645Corrected Total       199    19507.50Root MSE             7.64503    R-Square       0.40977Dependent Mean      51.85000    Adj R-Sq       0.40378Coeff Var           14.74451                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     18.11813    3.167133       5.72      <.0001math              1     0.663190    0.057872      11.46      <.0001female            1     -2.16840    1.086043      -2.00      0.0472

And here is OLS estimate for the second model.

The SYSLIN ProcedureOrdinary Least Squares EstimationModel                    WRITEDependent Variable       write                         Analysis of Variance                               Sum of        MeanSource                 DF     Squares      Square    F Value    Pr > FModel                   2    7856.321    3928.161      77.21    <.0001Error                 197    10022.55    50.87591Corrected Total       199    17878.88Root MSE             7.13273    R-Square       0.43942Dependent Mean      52.77500    Adj R-Sq       0.43373Coeff Var           13.51537                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     20.22837    2.713756       7.45      <.0001read              1     0.565887    0.049385      11.46      <.0001female            1     5.486894    1.014261       5.41      <.0001

Proc syslin with option sur also gives an estimate of the correlation between the errors of the two models. Here is the corresponding output.

The SYSLIN ProcedureSeemingly Unrelated Regression Estimation       Cross Model Covariance              SCIENCE       WRITESCIENCE        58.4464        7.8908WRITE           7.8908       50.8759      Cross Model Correlation              SCIENCE       WRITESCIENCE        1.00000       0.14471WRITE          0.14471       1.00000  Cross Model Inverse Correlation              SCIENCE       WRITESCIENCE        1.02139      -0.14780WRITE         -0.14780       1.02139   Cross Model Inverse Covariance              SCIENCE       WRITESCIENCE       0.017476      -.002710WRITE         -.002710      0.020076System Weighted MSE            0.9981Degrees of freedom                394System Weighted R-Square       0.3875

Finally, we have the seemingly unrelated regression estimation for our models. Note that both the estimates of the coefficients and their standard errors are different from the OLS model estimates shown above. 

The SYSLIN ProcedureSeemingly Unrelated Regression EstimationModel                  SCIENCEDependent Variable     science                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     20.13265    3.149485       6.39      <.0001math              1     0.625141    0.057528      10.87      <.0001female            1     -2.18934    1.086038      -2.02      0.0452Model                    WRITEDependent Variable       write                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     21.83439    2.698827       8.09      <.0001read              1     0.535484    0.049091      10.91      <.0001female            1     5.453748    1.014244       5.38      <.0001

Now that we have estimated our models let's test the predictor variables. The test for female combines information from both models. The tests for math and read are actually equivalent to the t-tests above except that the results are displayed as F-tests.

proc syslin data = "c:sasreghsb2" sur ;   science:  model science = math female ;   write:    model write = read female ;   female: stest science.female = write.female =0;   math: stest science.math = 0;   read: stest write.read = 0; run;

      Test Results for Variable FEAMLE  Num DF      Den DF    F Value    Pr > F       2         394      18.48    0.0001      Test Results for Variable MATH  Num DF      Den DF    F Value    Pr > F       1         394     118.31    0.0001      Test Results for Variable READ  Num DF      Den DF    F Value    Pr > F       1         394     119.21    0.0001

Now, let's estimate 3 models where we use the same predictors in each model as shown below.

     read  = female prog1 prog3      write = female prog1 prog3      math  = female prog1 prog3

Here variable prog1 and prog3 are dummy variables for the variable prog. Let's generate these variables before estimating our three models using proc syslin. 

data hsb2;   set "c:sasreghsb2";   prog1 = (prog = 1);   prog3 = (prog = 3); run;  proc syslin data = hsb2 sur;   model1: model read  = female prog1 prog3;   model2: model write = female prog1 prog3;   model3: model  math  = female prog1 prog3; run;

The OLS regression estimate of our three models are as follows. 

<some output omitted>

The SYSLIN ProcedureOrdinary Least Squares EstimationModel                   MODEL1Dependent Variable        read                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     56.82950    1.170562      48.55      <.0001female            1     -1.20858    1.327672      -0.91      0.3638prog1             1     -6.42937    1.665893      -3.86      0.0002prog3             1     -9.97687    1.606428      -6.21      <.0001Model                   MODEL2Dependent Variable       write                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     53.62162    1.042019      51.46      <.0001female            1     4.771211    1.181876       4.04      <.0001prog1             1     -4.83293    1.482956      -3.26      0.0013prog3             1     -9.43807    1.430021      -6.60      <.0001Model                   MODEL3Dependent Variable        math                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     57.10551    1.036890      55.07      <.0001female            1     -0.67377    1.176059      -0.57      0.5674prog1             1     -6.72394    1.475657      -4.56      <.0001prog3             1     -10.3217    1.422983      -7.25      <.0001

These regressions provide fine estimates of the coefficients and standard errors but these results assume the residuals of each analysis are completely independent of the other. Also, if we wish to test female, we would have to do it three times and would not be able to combine the information from all three tests into a single overall test.

Now let's see the output of the estimate using seemingly unrelated regression.

The SYSLIN ProcedureSeemingly Unrelated Regression Estimation

Model                   MODEL1Dependent Variable       read

                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     56.82950    1.170562      48.55      <.0001female            1     -1.20858    1.327672      -0.91      0.3638prog1             1     -6.42937    1.665893      -3.86      0.0002prog3             1     -9.97687    1.606428      -6.21      <.0001Model                   MODEL2Dependent Variable       write                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     53.62162    1.042019      51.46      <.0001female            1     4.771211    1.181876       4.04      <.0001prog1             1     -4.83293    1.482956      -3.26      0.0013prog3             1     -9.43807    1.430021      -6.60      <.0001Model                   MODEL3Dependent Variable        math                        Parameter Estimates                       Parameter    StandardVariable         DF     Estimate       Error    t Value    Pr > |t|Intercept         1     57.10551    1.036890      55.07      <.0001female            1     -0.67377    1.176059      -0.57      0.5674prog1             1     -6.72394    1.475657      -4.56      <.0001prog3             1     -10.3217    1.422983      -7.25      <.0001

Note that the coefficients are identical in the OLS results and in the seemingly unrelated regression estimate, however the standard errors are different, only slightly, due to the correlation among the residuals in the multiple equations. 

In addition to getting more appropriate standard errors, we can test the effects of the predictors across the equations.  We can test the hypothesis that the coefficient for female is 0 for all three outcome variables, as shown below.

proc syslin data = hsb2 sur;   model1: model read  = female prog1 prog3;   model2: model write = female prog1 prog3;   model3: model  math  = female prog1 prog3;   feamle: stest model1.female = model2.female = model3.female = 0; run;

     Test Results for Variable FEAMLE  Num DF      Den DF    F Value    Pr > F       3         588      11.63    0.0001

We can also test the hypothesis that the coefficient for female is 0 for just read and math. 

proc syslin data = hsb2 sur;   model1: model read  = female prog1 prog3;   model2: model write = female prog1 prog3;   model3: model  math  = female prog1 prog3;   f1: stest model1.female = model3.female = 0; run;

       Test Results for Variable F1  Num DF      Den DF    F Value    Pr > F       2         588       0.42    0.6599

We can also test the hypothesis that the coefficients for prog1 and prog3 are 0 for all three outcome variables, as shown below.

proc syslin data = hsb2 sur;   model1: model read  = female prog1 prog3;   model2: model write = female prog1 prog3;   model3: model  math  = female prog1 prog3;   progs: stest model1.prog1 = model2.prog1 = model3.prog1 =0,                model1.prog3 = model2.prog3 = model3.prog3 =0 ; run;

     Test Results for Variable PROGS  Num DF      Den DF    F Value    Pr > F       6         588      11.83    0.0001

4.5.2 Multivariate Regression

Let's now use multivariate regression using proc reg  to look at the same analysis say that we saw in the proc syslin example above, estimating the following 3 models.

     read  = female prog1 prog3      write = female prog1 prog3      math  = female prog1 prog3

Below we use proc reg to predict read write and math from female prog1 and prog3. Note that the top part of the output is similar to the sureg output in that it gives an overall summary of the model for each outcome variable, however the results are somewhat different and the sureg uses a Chi-Square test for the overall fit of the model, and mvreg uses an F-test. The lower part of the output appears similar to the sureg output, however when you compare the standard errors you see that the results are not the same.  These standard errors correspond to the OLS standard errors, so these results below do not take into account the correlations among the residuals (as do the sureg results).

proc reg data =hsb2;   model read write math = female prog1 prog3 ; run;

The REG Procedure [Some output omitted] Dependent Variable: read                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1       56.82950        1.17056      48.55      <.0001female        1       -1.20858        1.32767      -0.91      0.3638prog1         1       -6.42937        1.66589      -3.86      0.0002prog3         1       -9.97687        1.60643      -6.21      <.0001Dependent Variable: write                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1       53.62162        1.04202      51.46      <.0001female        1        4.77121        1.18188       4.04      <.0001prog1         1       -4.83293        1.48296      -3.26      0.0013prog3         1       -9.43807        1.43002      -6.60      <.0001Dependent Variable: math                        Parameter Estimates                     Parameter       StandardVariable     DF       Estimate          Error    t Value    Pr > |t|Intercept     1       57.10551        1.03689      55.07      <.0001female        1       -0.67377        1.17606      -0.57      0.5674prog1         1       -6.72394        1.47566      -4.56      <.0001prog3         1      -10.32168        1.42298      -7.25      <.0001

Now, let's test female. Note, that female was statistically significant in only one of the three equations. Using the mtest statement after proc reg allows us to test female across all three equations simultaneously. And, guess what? It is significant. This is consistent with what we found using seemingly unrelated regression estimation.

    female: mtest female=0;   run;

Multivariate Test: female                Multivariate Statistics and Exact F Statistics                             S=1    M=0.5    N=96Statistic                        Value    F Value    Num DF    Den DF    Pr > FWilks' Lambda               0.84892448      11.51         3       194    <.0001Pillai's Trace              0.15107552      11.51         3       194    <.0001Hotelling-Lawley Trace      0.17796108      11.51         3       194    <.0001Roy's Greatest Root         0.17796108      11.51         3       194    <.0001

We can also test prog1 and prog3, both separately and combined. Remember these are multivariate tests.

  prog1: mtest prog1 = 0; run;

Multivariate Test: prog1                Multivariate Statistics and Exact F Statistics                             S=1    M=0.5    N=96Statistic                        Value    F Value    Num DF    Den DF    Pr > FWilks' Lambda               0.89429287       7.64         3       194    <.0001Pillai's Trace              0.10570713       7.64         3       194    <.0001Hotelling-Lawley Trace      0.11820192       7.64         3       194    <.0001Roy's Greatest Root         0.11820192       7.64         3       194    <.0001

  prog3: mtest prog3 = 0; run;

Multivariate Test: prog3                Multivariate Statistics and Exact F Statistics                             S=1    M=0.5    N=96Statistic                        Value    F Value    Num DF    Den DF    Pr > FWilks' Lambda               0.75267026      21.25         3       194    <.0001Pillai's Trace              0.24732974      21.25         3       194    <.0001Hotelling-Lawley Trace      0.32860304      21.25         3       194    <.0001Roy's Greatest Root         0.32860304      21.25         3       194    <.0001

   prog: mtest prog1 = prog3 =0; run; quit;

Multivariate Test: prog                 Multivariate Statistics and F Approximations                              S=2    M=0    N=96Statistic                        Value    F Value    Num DF    Den DF    Pr > FWilks' Lambda               0.73294667      10.87         6       388    <.0001Pillai's Trace              0.26859190      10.08         6       390    <.0001Hotelling-Lawley Trace      0.36225660      11.68         6     256.9    <.0001Roy's Greatest Root         0.35636617      23.16         3       195    <.0001NOTE: F Statistic for Roy's Greatest Root is an upper bound.NOTE: F Statistic for Wilks' Lambda is exact.

Proc syslin with sur option and proc reg both allow you to test multi-equation models while taking into account the fact that the equations are not independent.  The proc syslin  with sur option allows you to get estimates for each equation which adjust for the non-independence of the equations, and it allows you to estimate equations which don't necessarily have the same predictors. By contrast, proc reg is restricted to equations that have the same set of predictors, and the estimates it provides for the individual equations are the same as the OLS estimates.  However, proc reg  allows you to perform more traditional multivariate tests of predictors.

4.6 Summary 

This chapter has covered a variety of topics that go beyond ordinary least squares regression, but there still remain a variety of topics we wish we could have covered, including the analysis of survey data, dealing with missing data, panel data analysis, and more. And, for the topics we did cover, we wish we could have gone into even more detail. One of our main goals for this chapter was to help you be aware of some of the techniques that are available in SAS for analyzing data that do not fit the assumptions of OLS regression and some of the remedies that are possible. If you are a member of the UCLA research community, and you have further questions, we invite you to use our consulting services to discuss issues specific to your data analysis.