Sample 25023: %VARMOD macro to fit a normal model accommodating unequal variances
 |  |  |  |  |  |
%VARMOD macro to fit a normal model accommodating unequal variances
| Contents: | Purpose / Requirements / Usage / Details / Limitations / Missing Values / See Also / References |
| NOTE: | Beginning in SAS 8.1, PROC QLIM in SAS/ETS software makes this macro obsolete. Use the HETERO statement in PROC QLIM to fit the model as described and illustrated in this Sample. |
- PURPOSE:
- The %VARMOD macro fits a model that accommodates heteroscedasticity (unequal variances) by iteratively fitting a normal regression model for the mean response and a log-linear regression model for the variance. Different linear functions can be used for the mean and log-variance models. The macro iteratively fits the two models until convergence of the likelihood function is obtained. Parameter estimates of both the mean and variance models are presented as well as a likelihood ratio test of the improvement of the fitted model over an equal variances model. Output data sets can be created containing parameter estimates, diagnostics and predicted values.
- REQUIREMENTS:
- The %VARMOD macro requires SAS/STAT Software, Version 7 or later. PROC QLIM requires SAS/ETS Software, Version 8.1 or later. The example presented uses syntax available beginning in SAS 9.0.
- USAGE:
-
Follow the instructions in the Downloads tab of this
sample to save the %VARMOD macro definition. Replace the text within quotes in the following statement with the location of the %VARMOD macro definition file on your system. In your SAS program or in the SAS editor window, specify this statement to define the %VARMOD macro and make it available for use:
%inc "<location of your file containing the VARMOD macro>";
Following this statement, you may call the %VARMOD macro. See the Results tab for an example.
The RESPONSE= option must be specified when using the macro to identify the response variable in the DATA= data set. Generally, you will also use the MMODEL= and VMODEL= options to specify the mean and variance models, respectively. You can use the same syntax used on the right-hand side of the MODEL statement in the GLM procedure. If either is omitted, an intercept-only model is fit. If either model involves CLASS variables, list the variables in the MCLASS= and/or VCLASS= options as needed. Several output data sets may be requested including data sets containing parameter estimates, diagnostics and predicted values. Predicted values can be obtained for a separate data set using the SCORE and OUTSCORE= options. There are a few additional options that may be needed occasionally.
Following is a list of available options. Note that the RESPONSE= option is required.
Options to identify input data set, response variable, and data set to score
- data=
- Specify the data set to be used. If omitted, the last-created data set is used.
- response=
- [REQUIRED] The response variable for the mean model.
- score=
- Specify an optional data set of observations for which predicted means and variances will be computed. These observations do not affect the model fit. The response variable, if present, is copied to the variable _OBSRESP. All explanatory variable values must be nonmissing for predicted values to be computed. The OUTSCORE= option should also be specified when SCORE= is specified.
Options to specify mean and variance models
- mclass=
- Specify all CLASS variables (if any) used in the mean model.
- mmodel=
- Specify all terms in the mean model using PROC GLM MODEL statement syntax.
- vclass=
- Specify all CLASS variables (if any) used in the variance model.
- vmodel=
- Specify all terms in the variance model using PROC GLM MODEL statement syntax.
- mnoint=
- If non-blank, fit a no-intercept mean model.
Options for output data sets
- outmdiag=
- Specify name of output data set of predicted values of the mean and diagnostics.
- outvdiag=
- Specify name of output data set of predicted values of the variance and diagnostics.
- outpred=
- Specify name of output data set containing the DATA= data set plus a variable _PREDM giving the predicted mean values and a variable _PREDV giving the predicted variances.
- outscore=
- Specify name of output data set containing the SCORE= data set plus a variable _PREDM giving the predicted mean values and a variable _PREDV giving the predicted variances. SCORE= should be specified when OUTSCORE= is specified.
- outmparm=
- Specify name of output data set containing parameter estimates of the mean model.
- outvparm=
- Specify name of output data set containing parameter estimates of the variance model.
- outlr=
- Specify name of output data set containing the likelihood ratio test for equal variances.
Options to control printing
- noprint=
- If non-blank, suppress all printed output. Useful when you just want to create output data sets.
- noeqvar=
- If non-blank, suppress printing equal variances mean model results.
- noiter=
- If non-blank, suppress printing iteration history.
- nolr=
- If non-blank, suppress printing likelihood ratio test for equal variances.
- noparms=
- If non-blank, suppress printing parameters of the mean and variance models.
- mgof=
- If non-blank, prints goodness of fit table for final mean model.
- vgof=
- If non-blank, prints goodness of fit table for final variance model.
Options to control maximum likelihood process
- maxiter=
- Maximum number of iterations allowed. An iteration is one fitting of both models. Default is 20.
- converge=
- Convergence is declared when the maximum change in loglikelihood is less than this. Default is 0.001.
- DETAILS:
-
This macro generalizes and enhances the %SEMIMOD macro that originally
appeared in the "Modeling Variance Homogeneity" example in the
PROC GENMOD documentation.
Note that the transformation presented by Box and Cox (1964) is also a possible solution to heteroscedasticity. This transformation can be done in PROC TRANSREG or PROC QLIM beginning in SAS 8.2.
For Taguchi designs in which there are several measurements taken at each setting of the control factors, this macro can be particularly useful since both the mean and the variance can be modeled. You can find which control factors affect the mean, which affect the variance and which affect both. You can then determine settings of the control factors that provide optimal mean response and variance.
The model fit by the %VARMOD macro is:
- y = Xβ + e , and
- V(y) = exp(ZL) , or equivalently: log(V(y)) = ZL ,
where y is the vector of responses, X is the model matrix for the mean model, β is the vector of parameter estimates for the mean model, e is a normal random variable with mean zero and variance V(y), Z is the model matrix for the variance model, and L is the vector of parameter estimates for the variance model. Note that the explanatory variables represented by the columns of X and Z are often common to the two models, but need not be.
Parameters of the models are estimated using an iterative maximum likelihood method which starts with an equal variances model. A likelihood ratio test of equal variances is provided using the likelihood values of the equal variances model and final fitted model. Small p-values indicate that at least one of the parameters in the variance model is significant and that the variances are not equal.
For each variance model parameter, a variance ratio is also computed. For the intercept, this is an estimate of the overall variance given that all other variance parameters are zero. For main effects, the variance ratio gives the change in variance comparing the current level to the reference level (for CLASS variables) or for a unit increase in the variable (for continuous variables). A zero parameter estimate corresponds to a variance ratio of one. Hence, the test of the parameter estimate is also a test of the variance ratio. The variance ratio for a change in x units of the variable can be computed as the variance ratio raised to the power x.
%VARMOD does not insist that the likelihood strictly increase. It simply looks for the change in successive iterations to be less than the CONVERGE= option value.
- LIMITATIONS:
-
Limited error checking is done. Be sure that:
- the DATA= data set exists and that the variables specified in the RESPONSE=, MCLASS=, MMODEL=, VCLASS=, and VMODEL= options exist on that data set
- the SCORE= data set exists if specified
- option names are specified as shown above
- the variables in your input (DATA=) data set do not have any of the following names: YVAR1 PRED XBETA STD HESSWGT UPPER LOWER RESRAW RESCHI RESDEV _SUM _WGT _ITER _DEVOLD _RSQUARE
- prior to invoking %VARMOD you do not have any data sets you want to keep with any of the following names: _TMP _WORK _OLD _ITER _MNFIT _MNDIAG _MEANEST _VARFIT _VARDIAG _VARESTS _LREQV _SCORE
- MISSING VALUES:
- Observations with missing values for any of the explanatory variables in the mean or variance model are omitted from the analysis. Observations with missing response, but which are nonmissing on all explanatory variables, are ignored in the model fitting phase, but predicted values are generated for these observations in the OUTMDIAG=, OUTVDIAG=, and OUTPRED= data sets. However, it is generally more convenient to use the SCORE= and OUTSCORE= options to get predicted values for new observations.
- SEE ALSO:
- Documentation for the QLIM procedure in the SAS/ETS User's Guide.
- REFERENCES:
-
Aitkin, M. (1987), "Modelling Variance Heterogeneity in Normal
Regression Using GLIM," Applied Statistics, 36(3), 332-339.
NOTE: Differences between -2*loglikelihood in %VARMOD and Aitkin's "deviance" is mostly due to his addition of log(2π) in the likelihood which is unnecessary since it is a constant.Box, G.E.P. and Cox, D.R. (1964), "An Analysis of Transformations," Journal of the Royal Statistics Society B, 26, 211-252.
Cook, R.D. and Weisberg, S. (1983), "Diagnostics for Heteroscedasticity in Regression," Biometrika, 70, 1-10.
In each of the Examples which follow, the analysis is done using both the %VARMOD macro and PROC QLIM. In order obtain the intercept for the variance model and the variance ratios, the following macro is defined. We'll reuse this macro throughout these examples to complete the QLIM analyses. Note that the macro needs the ParameterEstimates table from PROC QLIM to perform its computations. This requires that you specify an ODS OUTPUT statement to output this table when running PROC QLIM.
%macro qlimvarratios(inest=_last_);
* Compute the variance ratios and the variance model intercept ;
data _pe;
set &inest;
varratio=exp(estimate);
if substr(parameter,1,2) ne '_H' then varratio=.;
label varratio="Variance Ratio";
if parameter="_Sigma" then do;
output;
parameter="_H.Int";
* Var(intercept) determined by delta method ;
stderr=sqrt( 4*stderr**2/estimate**2 );
* VARMOD intercept = log(sigma**2) from QLIM ;
estimate=log(estimate**2);
if stderr>0 then do;
tvalue=estimate/stderr;
probt=2*probnorm(-abs(tvalue));
end;
varratio=exp(estimate);
end;
output;
run;
proc print data=_pe label noobs;
run;
%mend;
- EXAMPLE 1: Hypothetical semiconductor production
-
The data are generated from two hypothetical semiconductor
processes. The output, Y, of each process is affected by a
continuous process variable, X. Note that Process B has larger
output variability.
- Process A: Y = 1X + N(0,4)
- Process B: Y = 3X + N(0,25)
The variance ratio is 4/25 = 0.16. A single model for the mean is:
- E(Y) = 1*p1*X + 3*p2*X ,
where p1=1 if Process A and 0 otherwise; p2=1 if Process B and 0 otherwise. So, the mean model is simply the interaction term P*X. The variance model is:
- log(V(Y)) = log(4)*p1 + log(25)*p2 = 1.3863*p1 + 3.2189*p2 .
In a model with intercept made identifiable by setting the last parameter estimate to zero (as GENMOD and GLM do when a CLASS statement is specified), the variance model becomes:
- log(V(Y)) = log(25) + [log(4)-log(25)]*p1 + 0*p2 = 3.2189 + -1.8326*p1
This is the form of the variance model fit by %VARMOD when CLASS variable P is specified as the variance model. The variance ratio is then just exp(-1.8326) = 0.16. In PROC QLIM when the LINK=EXP and NOCONST options are specified in the HETERO statement, the variance model is of the form V(Y) = _Sigma2 * exp(ZL) and Z does not contain a column of ones so that L does not contain an intercept parameter. However, the QLIM variance model can be rewritten as:
- log(V(Y)) = log(_Sigma2) + ZL
which shows that log(_Sigma2) is the intercept in the variance model. QLIM provides estimates of _Sigma and the parameters in L (identified by the _H. prefix). For this example, Z need only contain the p1 column. Using the CLASS variable P in the HETERO statement would overparameterize the variance model.
The following statements generate the data and fit the mean and variance models as discussed above:
data semi; drop j; do p = 'a','b'; do x = 10 to 15; if p = 'a' then do; do j = 1 to 100; y = x + 2*rannor(13020238); output; end; end; else if p = 'b' then do; do j = 1 to 100; y = 3*x + 5*rannor(45678421); output; end; end; end; end; run; /* Define the VARMOD macro */ %inc "<location of your file containing the VARMOD macro>"; %varmod(data=semi, response=y, mclass=p, mmodel=p*x, vclass=p, vmodel=p) * Do the same analysis with PROC QLIM using syntax available in SAS 9.0 * or later *======================================================================; data semi2; set semi; p1=(p='a'); run; ods output parameterestimates=pe; proc qlim data=semi2; class p; model y = p*x; hetero y ~ p1 / link=exp noconst; run; %qlimvarratios; - RESULTS:
-
The iteration history is displayed in the log:
VARMOD Iteration 0: Fitting mean model VARMOD Iteration 0: Fitting variance model VARMOD Iteration 0: -2*LogLikelihood = 4339.4896598 VARMOD Iteration 1: Fitting mean model VARMOD Iteration 1: Fitting variance model VARMOD Iteration 1: -2*LogLikelihood = 3867.5374649 VARMOD Iteration 2: Fitting mean model VARMOD Iteration 2: Fitting variance model VARMOD Iteration 2: -2*LogLikelihood = 3867.1559993 VARMOD Iteration 3: Fitting mean model VARMOD Iteration 3: Fitting variance model VARMOD Iteration 3: -2*LogLikelihood = 3867.1559992 VARMOD: Convergence attained.
Following are the results of the analysis. The likelihood ratio test shows that the fitted model, which allows for variance differences, is significantly better than an equal variances model (p<.0001). Note that the likelihood ratio chi-square value is the difference in -2*log(likelihood) values at the initial (equal variances) and final iterations. To obtain this test using PROC QLIM, run QLIM twice: once with the HETERO statement and once without (the equal variances model). Then compute twice the difference in the log likelihoods for the two models.
The estimated parameters of the final mean and variance models closely match the true parameters that generated the data. In the QLIM results, the effects preceding _Sigma are the mean model parameters and the effects following _Sigma are the variance model parameters. The variance ratio of 23.79 for the variance model intercept represents the estimated variance when all variance model parameters are zero, which corresponds to Process B whose true variance is 25. The variance ratio for P represents the estimated ratio of process variances (A to B). The estimated value of 0.15 also closely matches the true variance ratio of 0.16. Since this indicates that the process B variance is larger than the process A variance, you may prefer to use the B to A ratio which is simply the reciprocal, 1/0.1506 = 6.64.
Equal variances model The GENMOD Procedure Model Information Data Set WORK._WORK Distribution Normal Link Function Identity Dependent Variable y Scale Weight Variable _wgt Number of Observations Read 1200 Number of Observations Used 1200 Sum of Weights 1200 Class Level Information Class Levels Values p 2 a b Criteria For Assessing Goodness Of Fit Criterion DF Value Value/DF Deviance 1197 16421.0318 13.7185 Scaled Deviance 1197 16421.0318 13.7185 Pearson Chi-Square 1197 16421.0318 13.7185 Scaled Pearson X2 1197 16421.0318 13.7185 Log Likelihood -9313.2421 Algorithm converged. Analysis Of Parameter Estimates Standard Wald 95% Chi- Parameter DF Estimate Error Confidence Limits Square Pr > ChiSq Intercept 1 -1.0176 0.2133 -1.4355 -0.5996 22.77 <.0001 x*p a 1 1.0880 0.0171 1.0546 1.1215 4068.87 <.0001 x*p b 1 3.0920 0.0171 3.0585 3.1254 32859.0 <.0001 Scale 0 1.0000 0.0000 1.0000 1.0000 NOTE: The scale parameter was held fixed. Equal variances model The GENMOD Procedure Lagrange Multiplier Statistics Parameter Chi-Square Pr > ChiSq Scale 4820.3264 <.0001 Iteration History -2*Log Iteration Likelihood 0 4339.4896597855 1 3867.537464924 2 3867.1559992645 3 3867.1559992075 Likelihood Ratio Test of Equal Variances LR Chi Square DF Pr>Chi 472.334 1 <.0001 Final Mean Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Parameter Level1 DF Estimate Error Limit Limit Square ChiSq Intercept 1 -0.6888 0.5321 -1.7318 0.3542 1.68 0.1955 x*p a 1 1.0622 0.0422 0.9794 1.1450 632.46 <.0001 x*p b 1 3.0662 0.0447 2.9786 3.1537 4710.98 <.0001 Scale 0 1.0000 0.0000 1.0000 1.0000 _ _ Final Variance Model Analysis of Parameter Estimates 95% Lower Standard Confidence Parameter Level1 DF Estimate Error Limit Intercept 1 3.1693 0.0577 3.0561 p a 1 -1.8933 0.0816 -2.0533 p b 0 0.0000 0.0000 0.0000 Scale 0 0.5000 0.0000 0.5000 95% Upper Confidence Chi- Pr > Variance Limit Square ChiSq Ratio 3.2824 3013.28 <.0001 23.7901 -1.7332 537.68 <.0001 0.1506 0.0000 . . . 0.5000 _ _ . The QLIM Procedure Summary Statistics of Continuous Responses N Obs N Obs Standard Lower Upper Lower Upper Variable Mean Error Type Bound Bound Bound Bound y 25.10754 13.651861 Regular Class Level Information Class Levels Values p 2 a b Model Fit Summary Number of Endogenous Variables 1 Endogenous Variable y Number of Observations 1200 Log Likelihood -3036 Maximum Absolute Gradient 0.0004995 Number of Iterations 12 AIC 6083 Schwarz Criterion 6108 Algorithm converged. Parameter Estimates Standard Approx Parameter Estimate Error t Value Pr > |t| Intercept -0.688780 0.532205 -1.29 0.1956 x*p a 1.062221 0.042243 25.15 <.0001 x*p b 3.066155 0.044677 68.63 <.0001 _Sigma 4.877509 0.140810 34.64 <.0001 _H.p1 -1.893279 0.081660 -23.19 <.0001 QLIM results Standard Approx Variance Parameter Level1 Estimate Error t Value Pr > |t| Ratio Intercept -0.688780 0.532205 -1.29 0.1956 . x*p a 1.062221 0.042243 25.15 <.0001 . x*p b 3.066155 0.044677 68.63 <.0001 . _Sigma 4.877509 0.140810 34.64 <.0001 . _H.Int 3.169269 0.057739 54.89 <.0001 23.7901 _H.p1 -1.893279 0.081660 -23.19 <.0001 0.1506 - EXAMPLE 2:
-
Reference: Byrne, D.M. and Taguchi, G. (1986), "The Taguchi approach to parameter design," ASQC 40th Anniversary
Quality Control Congress Transactions, Milwaukee, WI: American
Society of Quality Control, 168-177.
data bt; input SHRINK TIME MTEMP THICK HPRESS SPEED HTIME GATE; datalines; 2.2 -1 -1 -1 -1 -1 -1 -1 2.1 -1 -1 -1 -1 -1 -1 -1 2.3 -1 -1 -1 -1 -1 -1 -1 2.3 -1 -1 -1 -1 -1 -1 -1 0.3 -1 -1 -1 1 1 1 1 2.5 -1 -1 -1 1 1 1 1 2.7 -1 -1 -1 1 1 1 1 0.3 -1 -1 -1 1 1 1 1 0.5 -1 1 1 -1 -1 1 1 3.1 -1 1 1 -1 -1 1 1 0.4 -1 1 1 -1 -1 1 1 2.8 -1 1 1 -1 -1 1 1 2 -1 1 1 1 1 -1 -1 1.9 -1 1 1 1 1 -1 -1 1.8 -1 1 1 1 1 -1 -1 2 -1 1 1 1 1 -1 -1 3 1 -1 1 -1 1 -1 1 3.1 1 -1 1 -1 1 -1 1 3 1 -1 1 -1 1 -1 1 3 1 -1 1 -1 1 -1 1 2.1 1 -1 1 1 -1 1 -1 4.2 1 -1 1 1 -1 1 -1 1 1 -1 1 1 -1 1 -1 3.1 1 -1 1 1 -1 1 -1 4 1 1 -1 -1 1 1 -1 1.9 1 1 -1 -1 1 1 -1 4.6 1 1 -1 -1 1 1 -1 2.2 1 1 -1 -1 1 1 -1 2 1 1 -1 1 -1 -1 1 1.9 1 1 -1 1 -1 -1 1 1.9 1 1 -1 1 -1 -1 1 1.8 1 1 -1 1 -1 -1 1 ;Notice the affect on SHRINK variability due to HTIME shown by the side-by-side box plots produced by PROC BOXPLOT:proc sort; by htime; run; proc boxplot; plot shrink*htime/ cframe = vligb cboxes = dagr cboxfill = ywh; run;
Fit a model allowing variance heterogeneity due to HTIME.
%varmod(data=bt, response=shrink, mmodel=time mtemp thick hpress speed htime gate, vmodel=htime)The iteration history is displayed in the log:VARMOD Iteration 0: Fitting mean model VARMOD Iteration 0: Fitting variance model VARMOD Iteration 0: -2*LogLikelihood = 20.870848111 VARMOD Iteration 1: Fitting mean model VARMOD Iteration 1: Fitting variance model VARMOD Iteration 1: -2*LogLikelihood = -46.81346632 VARMOD Iteration 2: Fitting mean model VARMOD Iteration 2: Fitting variance model VARMOD Iteration 2: -2*LogLikelihood = -46.81346632 VARMOD: Convergence attained.
Following are selected tables from the results of the macro. The likelihood ratio test shows that the model allowing different variances at the HTIME levels is an improvement over the equal variances model. In the Final Variance Model table, the variance ratio for HTIME is 16.52 which is the change in variance for a one unit increase in HTIME. Note that the coded values of the low and high levels of HTIME differ by two (-1 to 1). The variance ratio for a change in x units is computed by raising the variance ratio for one unit to the power x. So, the variance ratio comparing the low and high levels of HTIME requires squaring the printed value, resulting in 272.94. Hence, the variance at high HTIME is about 273 times the variance at low HTIME. This can be seen by using the OUTPRED= option to create a data set of predicted values. See the predicted variance values in the _PREDV variable of this data set.Likelihood Ratio Test of Equal Variances LR Chi Square DF Pr>Chi 67.6843 1 <.0001 Final Mean Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Parameter DF Estimate Error Limit Limit Square ChiSq Intercept 1 2.2500 0.1486 1.9588 2.5412 229.38 <.0001 TIME 1 0.4250 0.1486 0.1338 0.7162 8.18 0.0042 MTEMP 1 -0.0750 0.1486 -0.3662 0.2162 0.25 0.6137 THICK 1 0.0625 0.1486 -0.2287 0.3537 0.18 0.6740 HPRESS 1 -0.2813 0.1486 -0.5724 0.0099 3.58 0.0583 SPEED 1 0.1438 0.1486 -0.1474 0.4349 0.94 0.3332 HTIME 1 -0.0187 0.1486 -0.3099 0.2724 0.02 0.8996 GATE 1 -0.2313 0.1486 -0.5224 0.0599 2.42 0.1196 Scale 0 1.0000 0.0000 1.0000 1.0000 _ _ Final Variance Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Variance Parameter DF Estimate Error Limit Limit Square ChiSq Ratio Intercept 1 -2.4629 0.2500 -2.9529 -1.9729 97.06 <.0001 0.0852 HTIME 1 2.8046 0.2500 2.3146 3.2946 125.85 <.0001 16.5209 Scale 0 0.5000 0.0000 0.5000 0.5000 _ _ .The following statements perform the same analysis using PROC QLIM which produces the same results.ods output parameterestimates=pe; proc qlim data=bt; model shrink = time mtemp thick hpress speed htime gate; hetero shrink ~ htime / link=exp noconst; run; %qlimvarratios; Standard Approx Variance Parameter Estimate Error t Value Pr > |t| Ratio Intercept 2.250000 0.148561 15.15 <.0001 . TIME 0.425000 0.148561 2.86 0.0042 . MTEMP -0.075000 0.148561 -0.50 0.6137 . THICK 0.062500 0.148561 0.42 0.6740 . HPRESS -0.281250 0.148561 -1.89 0.0583 . SPEED 0.143750 0.148561 0.97 0.3332 . HTIME -0.018750 0.148561 -0.13 0.8996 . GATE -0.231250 0.148561 -1.56 0.1196 . _Sigma 0.291866 0.036483 8.00 <.0001 . _H.Int -2.462921 0.250000 -9.85 <.0001 0.0852 _H.HTIME 2.804625 0.250000 11.22 <.0001 16.5209 - EXAMPLE 3:
-
The following shows how the sequence of models in Table 2 of Aitkin (1987) can be fit
using the %VARMOD macro and PROC QLIM. Results from %VARMOD and PROC QLIM are not shown.
Note that the results of Aitkin's final model (and the last fit in the below %VARMOD and
PROC QLIM sequences) is
shown in the Appendix of the paper. After the first run, the equal variances
output and the iteration history are suppressed (by the NOITER=NO and NOEQVAR=NO options).
Note that the -2*loglikelihood reported by %VARMOD differs from
Aitkin's "deviance" mostly due to his inclusion of log(2π) in the
likelihood. This is unnecessary since it is constant.
data treevol; input diam ht vol @@; y=vol**(1/3); htsq=ht*ht; diamsq=diam*diam; ht_diam=ht*diam; lht=log(ht); ldiam=log(diam); lvol=log(vol); cards; 8.3 70 10.3 8.6 65 10.3 8.8 63 10.2 10.5 72 16.4 10.7 81 18.8 10.8 83 19.7 11.0 66 15.6 11.0 75 18.2 11.1 80 22.6 11.2 75 19.9 11.3 79 24.2 11.4 76 21.0 11.4 76 21.4 11.7 69 21.3 12.0 75 19.1 12.9 74 22.2 12.9 85 33.8 13.3 86 27.4 13.7 71 25.7 13.8 64 24.9 14.0 78 34.5 14.2 80 31.7 14.5 74 36.3 16.0 72 38.3 16.3 77 42.6 17.3 81 55.4 17.5 82 55.7 17.9 80 58.3 18.0 80 51.5 18.0 80 51.0 20.6 87 77.0 ; %varmod(data=treevol, response=y, mmodel=ht diam, vmodel=ht diam) %varmod(data=treevol, response=y, mmodel=ht diam, vmodel=ht, noeqvar=no, noiter=no) %varmod(data=treevol, response=y, mmodel=ht diam) %varmod(data=treevol, response=y, mmodel=ht diam, vmodel=ht|diam htsq diamsq, noeqvar=no, noiter=no) %varmod(data=treevol, response=y, mmodel=ht diam, vmodel=ht diam htsq diamsq, noeqvar=no, noiter=no) %varmod(data=treevol, response=y, mmodel=ht diam, vmodel=ht diam diamsq, noeqvar=no, noiter=no) %varmod(data=treevol, response=y, mmodel=ht diam, vmodel=diam diamsq, noeqvar=no, noiter=no) %varmod(data=treevol, response=lvol, mmodel=lht ldiam, vmodel=ldiam ldiam*ldiam, noeqvar=no, noiter=no)The following perform the analyses using PROC QLIM. Note that several of these models do not converge properly, particularly the ones with the more complex variance models. The VARMOD macro seems more successful in these cases.ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; hetero y ~ ht diam / link=exp noconst; run; %qlimvarratios; ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; hetero y ~ ht / link=exp noconst; run; %qlimvarratios; ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; run; %qlimvarratios; ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; hetero y ~ ht|diam htsq diamsq / link=exp noconst; run; %qlimvarratios; ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; hetero y ~ ht diam htsq diamsq / link=exp noconst; run; %qlimvarratios; ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; hetero y ~ ht diam diamsq / link=exp noconst; run; %qlimvarratios; ods output parameterestimates=pe; proc qlim data=treevol; model y = ht diam; hetero y ~ diam diamsq / link=exp noconst; run; %qlimvarratios; proc qlim data=treevol; model lvol = lht ldiam; hetero lvol ~ ldiam ldiam*ldiam / link=exp noconst; run; %qlimvarratios;
- EXAMPLE 4:
-
This example is from Desrochers and Ewing (1984), "Leaf spring free height
analysis using Taguchi methods," Second Supplier Symposium on Taguchi
Methods, pp. 38-47.
The following statements recreate the experimental design in a data set and add the observed responses at each design point.
proc factex; /* Control factors design (resolution IV) - inner array */ factors temp heattime trantime holdtime; model res=4; size design=8; output out=in; run; /* Noise factors design - outer array */ factors oiltemp; output out=inout designrep=in; run; /* Replication */ factors rep /nlev=3; output out=inoutrep designrep=inout; run; quit; /* Response data */ data response; input height @@; cards; 7.56 7.62 7.44 7.18 7.18 7.25 7.5 7.56 7.5 7.5 7.56 7.5 7.94 8 7.88 7.32 7.44 7.44 7.78 7.78 7.81 7.5 7.25 7.12 7.56 7.81 7.69 7.81 7.5 7.59 7.59 7.56 7.75 7.63 7.75 7.56 7.69 8.09 8.06 7.56 7.69 7.62 8.15 8.18 7.88 7.88 7.88 7.44 ; /* Merge the full design with the response data */ data taguchi; merge inoutrep response; drop rep; run;
For a first model, all variables were used in both the mean and variance models. Tests of the variance model parameters indicated that TEMP and HOLDTIME do not significantly affect the variance though they do affect the mean. Increasing either variable increases the mean but does not significantly affect variability. So, the following model is fit.The unequal variances model is significant compared to the equal variances model. In this simpler model, TRANTIME seems to affect the variance but not the mean, so it can be removed from the mean model. Note that increasing TRANTIME allows you to decrease variability without significantly changing the mean.
%varmod(data=taguchi, response=height, mmodel=temp heattime trantime holdtime, vmodel=heattime trantime)Likelihood Ratio Test of Equal Variances LR Chi Square DF Pr>Chi 16.7662 2 0.0002 Final Mean Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Parameter DF Estimate Error Limit Limit Square ChiSq Intercept 1 7.6306 0.0293 7.5732 7.6880 67921.7 <.0001 temp 1 0.1067 0.0195 0.0684 0.1449 29.88 <.0001 heattime 1 0.0809 0.0281 0.0259 0.1359 8.30 0.0040 trantime 1 0.0277 0.0204 -0.0122 0.0675 1.85 0.1742 holdtime 1 0.0468 0.0195 0.0086 0.0851 5.75 0.0164 Scale 0 1.0000 0.0000 1.0000 1.0000 _ _ Final Variance Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Variance Parameter DF Estimate Error Limit Limit Square ChiSq Ratio Intercept 1 -3.5963 0.2041 -3.9964 -3.1963 310.41 <.0001 0.02742 heattime 1 0.9729 0.2041 0.5729 1.3730 22.72 <.0001 2.64573 trantime 1 -0.4340 0.2041 -0.8340 -0.0339 4.52 0.0335 0.64794 Scale 0 0.5000 0.0000 0.5000 0.5000 _ _ .Performing the same analysis using PROC QLIM:ods output parameterestimates=pe; proc qlim data=taguchi; model height = temp heattime trantime holdtime; hetero height ~ heattime trantime / link=exp noconst; run; %qlimvarratios; Standard Approx Variance Parameter Estimate Error t Value Pr > |t| Ratio Intercept 7.630614 0.029431 259.27 <.0001 . temp 0.106667 0.019583 5.45 <.0001 . heattime 0.080889 0.028283 2.86 0.0042 . trantime 0.027656 0.020426 1.35 0.1757 . holdtime 0.046809 0.019583 2.39 0.0168 . _Sigma 0.165602 0.016902 9.80 <.0001 . _H.Int -3.596338 0.204124 -17.62 <.0001 0.02742 _H.heattime 0.972938 0.229888 4.23 <.0001 2.64571 _H.trantime -0.433976 0.230527 -1.88 0.0598 0.64793The following statements create a set of observations over the range of HEATTIMEs with HOLDTIME, TRANTIME, and TEMP set at their maximums (explained below). Predicted means and variances will be obtained for these points using the final model.data heat; holdtime=1; trantime=1; temp=1; do heattime=-1 to 1 by .1; output; end; run;Since TRANTIME was found to not significantly affect the mean, the following mean model removes it. An output data set of predicted values is created using the OUTPRED= option. Predicted values for the additional observations in the HEAT data set are also output to a data set for use later using the SCORE= and OUTSCORE= options.%varmod(data=taguchi, response=height, score=heat, outscore=heatpred, outpred=pred, mmodel=temp heattime holdtime, vmodel=heattime trantime) Final Mean Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Parameter DF Estimate Error Limit Limit Square ChiSq Intercept 1 7.6418 0.0281 7.5867 7.6969 73993.7 <.0001 temp 1 0.1070 0.0202 0.0674 0.1466 28.04 <.0001 heattime 1 0.0810 0.0281 0.0260 0.1361 8.32 0.0039 holdtime 1 0.0471 0.0202 0.0075 0.0867 5.44 0.0196 Scale 0 1.0000 0.0000 1.0000 1.0000 _ _ Final Variance Model Analysis of Parameter Estimates 95% Lower 95% Upper Standard Confidence Confidence Chi- Pr > Variance Parameter DF Estimate Error Limit Limit Square ChiSq Ratio Intercept 1 -3.5598 0.2041 -3.9599 -3.1597 304.13 <.0001 0.02844 heattime 1 0.9199 0.2041 0.5198 1.3199 20.31 <.0001 2.50892 trantime 1 -0.4248 0.2041 -0.8249 -0.0248 4.33 0.0374 0.65388 Scale 0 0.5000 0.0000 0.5000 0.5000 _ _ .Performing the same analysis using PROC QLIM:ods output parameterestimates=pe; proc qlim data=taguchi; model height = temp heattime holdtime; hetero height ~ heattime trantime / link=exp noconst; run; %qlimvarratios; Standard Approx Variance Parameter Estimate Error t Value Pr > |t| Ratio Intercept 7.641807 0.028241 270.59 <.0001 . temp 0.106971 0.020276 5.28 <.0001 . heattime 0.081023 0.028318 2.86 0.0042 . holdtime 0.047130 0.020276 2.32 0.0201 . _Sigma 0.168655 0.017213 9.80 <.0001 . _H.Int -3.559796 0.204124 -17.44 <.0001 0.02844 _H.heattime 0.919894 0.231265 3.98 <.0001 2.50902 _H.trantime -0.424903 0.239876 -1.77 0.0765 0.65383Print the original data set with predicted means (_PREDM) and variances (_PREDV).proc print data=pred noobs; var temp heattime trantime holdtime height _predm _predv; run; temp heattime trantime holdtime height _predm _predv -1 -1 -1 -1 7.56 7.4066759 0.0173388 -1 -1 -1 -1 7.62 7.4066759 0.0173388 -1 -1 -1 -1 7.44 7.4066759 0.0173388 -1 -1 -1 -1 7.18 7.4066759 0.0173388 -1 -1 -1 -1 7.18 7.4066759 0.0173388 -1 -1 -1 -1 7.25 7.4066759 0.0173388 -1 -1 1 1 7.50 7.500932 0.0074133 -1 -1 1 1 7.56 7.500932 0.0074133 -1 -1 1 1 7.50 7.500932 0.0074133 -1 -1 1 1 7.50 7.500932 0.0074133 -1 -1 1 1 7.56 7.500932 0.0074133 -1 -1 1 1 7.50 7.500932 0.0074133 -1 1 -1 1 7.94 7.6629901 0.109142 -1 1 -1 1 8.00 7.6629901 0.109142 -1 1 -1 1 7.88 7.6629901 0.109142 -1 1 -1 1 7.32 7.6629901 0.109142 -1 1 -1 1 7.44 7.6629901 0.109142 -1 1 -1 1 7.44 7.6629901 0.109142 -1 1 1 -1 7.78 7.568734 0.0466641 -1 1 1 -1 7.78 7.568734 0.0466641 -1 1 1 -1 7.81 7.568734 0.0466641 -1 1 1 -1 7.50 7.568734 0.0466641 -1 1 1 -1 7.25 7.568734 0.0466641 -1 1 1 -1 7.12 7.568734 0.0466641 1 -1 -1 1 7.56 7.7148699 0.0173388 1 -1 -1 1 7.81 7.7148699 0.0173388 1 -1 -1 1 7.69 7.7148699 0.0173388 1 -1 -1 1 7.81 7.7148699 0.0173388 1 -1 -1 1 7.50 7.7148699 0.0173388 1 -1 -1 1 7.59 7.7148699 0.0173388 1 -1 1 -1 7.59 7.6206137 0.0074133 1 -1 1 -1 7.56 7.6206137 0.0074133 1 -1 1 -1 7.75 7.6206137 0.0074133 1 -1 1 -1 7.63 7.6206137 0.0074133 1 -1 1 -1 7.75 7.6206137 0.0074133 1 -1 1 -1 7.56 7.6206137 0.0074133 1 1 -1 -1 7.69 7.7826718 0.109142 1 1 -1 -1 8.09 7.7826718 0.109142 1 1 -1 -1 8.06 7.7826718 0.109142 1 1 -1 -1 7.56 7.7826718 0.109142 1 1 -1 -1 7.69 7.7826718 0.109142 1 1 -1 -1 7.62 7.7826718 0.109142 1 1 1 1 8.15 7.876928 0.0466641 1 1 1 1 8.18 7.876928 0.0466641 1 1 1 1 7.88 7.876928 0.0466641 1 1 1 1 7.88 7.876928 0.0466641 1 1 1 1 7.88 7.876928 0.0466641 1 1 1 1 7.44 7.876928 0.0466641
Using the predicted values from the HEAT data set, these statements plot the predicted mean and variance over the range of HEATTIME. TEMP and HOLDTIME are fixed at their maximum values since their positive parameter estimates in the mean model imply this maximizes the response. TRANTIME is also fixed at its maximum value since its negative parameter in the log-variance model implies this minimizes variance. HEATTIME is the only common variable to the mean and variance models and increasing it increases the response but also increases the variance.The desired response (HEIGHT) is 8 inches. If the maximum tolerable variation is +/- 0.50 inch, then this corresponds to a standard deviation of about 0.5/3 = 0.1667 or a variance of 0.0278. Using the estimated mean and variance models, a variance of 0.0278 corresponds to a HEATTIME of about 0.436 resulting in a height of approximately 7.831 inches -- about 0.17 inch below target. A larger tolerated variance would allow you to get closer to the target value of 8 inches. However, if it is possible to increase either or both of TEMP and HOLDTIME above their current maximums, this would increase HEIGHT without increasing the variance (assuming that the same model holds outside the current variable ranges).
proc gplot data=heatpred; plot _predm*heattime="M" / haxis=-1 to -.2 by .2,0,.2 to 1 by .2 href=0.436 vref=7.831; plot2 _predv*heattime="V" / vref=0.0278 ; label _predm="MEAN" _predv="VARIANCE"; run;
Right-click on the link below and select Save to save
the %VARMOD macro definition
to a file. It is recommended that you name the file
varmod.sas.
| Type: | Sample |
| Topic: | SAS Reference ==> Procedures ==> QLIM Analytics ==> Econometrics Analytics ==> Categorical Data Analysis SAS Reference ==> Procedures ==> REG Analytics ==> Multivariate Analysis SAS Reference ==> Procedures ==> GENMOD SAS Reference ==> Procedures ==> GLM Analytics ==> Longitudinal Analysis Analytics ==> Regression Analytics ==> Analysis of Variance |
| Date Modified: | 2007-08-14 03:03:09 |
| Date Created: | 2005-01-13 15:03:57 |
Operating System and Release Information
| Product Family | Product | Host | SAS Release | |
| Starting | Ending | |||
| SAS System | SAS/STAT | All | 8 TS M0 | n/a |
| SAS System | SAS/ETS | All | 8 TS M0 | n/a |
