![]() | ![]() | ![]() | ![]() | ![]() |
| 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. |
%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
Options to specify mean and variance models
Options for output data sets
Options to control printing
Options to control maximum likelihood process
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:
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.
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;
The variance ratio is 4/25 = 0.16. A single model for the mean is:
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:
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:
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:
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;
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
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
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;
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.64793
The 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.65383
Print 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.0466641Using 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 |
| 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 |





