**********************************************************; * One-way ANOVA and nonparametric `one-way layouts': * (i) Variability of treatment means and the one-way ANOVA test * (ii) One-sample tests and the normality of residuals * (iii) The Kruskal-Wallis test: A nonparametric one-way layout * * Assume that we have measurements for K different `treatments' or * `treatment groups', with multiple observations from each * treatment group. * * The word `treatment' is used even when `treatment' makes no sense * literally, such as products from different vendors or * measurements made on different days. * * Do the means of the treatment groups vary significantly, allowing * for random variation within treatment groups? * * Data that is arranged in this way is called a `one-way layout'. * * The traditional statistical model for a one-way layout is * * (*) X_ij = mu_i + e_ij where * * e_ij are N(0,sigma^2) = sigma*N(0,1) * * Here X_ij is the j-th observation for the i-th treatment. * The model (*) means that the data from different treatment groups * have means mu_i that depend on the treatment. There are also * normally-distributed ``errors'' e_ij that have the same variance * across treatments. * * Given (*), one can test * * H_0: All treatment groups have the same mean Vs. * * H_1: Means of treatment groups are not the same * * using a classical `one-way ANOVA' test. If K=2, this is * equivalent to the classical (equal-variance) Student-t test. * * If we use a one-way ANOVA to test H_0, we should check that * * (i) the errors e_ij in (*) are normally distributed. We can do * this by checking that the RESIDUALS of the model (*) --- * that is, * * R_ij = X_ij - Xbar_i * * where Xbar_i is the mean of the observations in the i-th * treatment group --- are consistent with a normal distribution. * The residuals are estimates of the errors e_ij in (*). * * (ii) the variance of the errors does not vary with the treatment * group. This we can also check from the residuals. * * If the errors are not normal or the error variances depend on the * treatment group, then the P-value from the one-way ANOVA test * may not be valid for small samples. The one-way ANOVA test is * reasonably robust for large samples (for example, 10 or more * observations per treatment group) if treatment-group variances * are the same. In that case, (i) is less of a concern. * * We can check (ii) by plotting the data by treatment group, which * will also give us an idea of what the data looks like. We will * check (i) by applying a statistical test for normality to the * residuals and also by plotting the residuals on ``normal paper.'' * * A NONPARAMETRIC TEST: * * The model (*) can be generalized with `sigma*N(0,1)' is replaced * by an arbitrary error distribution, which is assumed to be the * same for all treatment groups. In that case, one can use the * Kruskal-Wallis (KW) test, which is a generalization of the * Wilcoxon rank-sum test for K=2. The KW test is almost as * powerful as the one-way ANOVA test for normally-distributed * data, and can be much more powerful if, for example, the * observations X_ij are exponentially distributed. * * The KW test is based on assigning ranks to all the data from all * treatments and then doing (essentially) a one-way ANOVA test on * the ranks. The KW test allows for the possibilility that the * distribution of treatment errors is not normal, but still * requires that the error distribution (and in particular, its * variance) be the same across treatments. * * AN EXAMPLE OF A ONE-WAY LAYOUT: * * In the example that follows, we assume that we have measured * the hardness of alloy samples that resulted from four different * manufacturing processes. Thus, in this case, `treatment' really * does mean `treatment'. * * In the data step below, we read data word-by-word and `retain' * the words A or B or C or D as treatment names. * * Note the `do end' pair of commands in the `else' statement. This * allows us to include more than one command in the `else' block. * The use of `output' within the do-end pair means that a record * is created ONLY IN that case. This means that we will not get * records with missing values of hardness when the treat=zz * variable changes. * * The final command `drop zz' in the data step tells SAS to skip * the `zz' column in the output, which also saves space. * **********************************************************; title 'Hardness of Alloys from four Different Processes -- YOURNAME'; options ls=75 ps=60 nocenter pageno=1; data alloys; retain treat; input zz$ @@; if zz='A' or zz='B' or zz='C' or zz='D' then treat=zz; else do; hardness=input(zz,12.0); output; end; drop zz; datalines; A 166 194 211 213 238 B 169 175 185 195 202 218 C 205 227 234 270 D 102 145 148 195 202 run; **********************************************************; * We can get an idea of the distribution of the values over the ; * four treatments from a `subgroup' horizontal bar chart, ; * and of relative variances from a `group' bar chart: ; * * SAS takes an entire page of output for each chart: We can improve * the appearance of the charts by temporarily decreasing the * output page size. * * The two bar charts could be done in one `proc chart': The purpose * of two `proc chart's is so that we can use different title lines. **********************************************************; options ps=40; * Contract page height for nicer plots; proc chart data=alloys; title2 'Horizontal Bar Charts for alloy hardnesses'; title3 "Bar chart by `subgroup'"; title4 '(Treatments are A,B,C,D)'; hbar hardness / subgroup=treat; run; proc chart data=alloys; title3 "Bar chart by `group'"; title4 '(Treatments are A,B,C,D)'; title5 'This output gives a clearer picture of the'; title6 ' relative treatment-group variances.'; hbar hardness / group=treat; run; options ps=60; * Restore page height to its usual value; **********************************************************; * Now we carry out the statistical analysis for a one-way ANOVA. * * `Proc glm' is a SAS routine that carries out a wide variety of * linear models procedures. * * `Means treat' tells SAS to display a table with treatment-group * means and standard deviations, similar to what we would get in * `proc means' output. This can also be used to visually check * variation in treatment-group standard deviations. * * The line `output p=means r=resid' tells SAS to add two * additional columns to the current data set: * * (i) the treatment means (`p' stands for `predicted'), to have * the variable name `means' and * * (ii) the residuals (`r' for `residual'), to have the variable * name `resid'. * * In general, a `residual' is an observed value minus the * predicted value, here minus the treatment mean. We will use * two extra variables later to see how well the data fits the * model assumptions. ***********************************************************; proc glm data=alloys; title2 'ONE-WAY ANOVA ANALYSIS'; title3 'Also display treatment means and save residuals'; class treat; model hardness=treat; output p=means r=resid; means treat; run; proc print noobs; * Noobs says to skip the OBS column; title2 'THE CURRENT DATASET AS SAS SEES IT'; title3 'Note the two extra columns for means and residuals'; run; **********************************************************; * In general, we want to test whether the residual variances vary * with the treatment group in any way, and whether the residuals * appear to be normally distributed. * * We can skip the first test in this case, since we have already * plotted the observations within each treatment group. In more * complex examples, we might need several different plots to * check whether residual variances depend on variables such * as which `treatment group'. * * We test the residuals as a group to see if they appear to be * normally distributed. * * `Proc univariate' carries out a wide variety of tests and * procedures on a single variable. (That is, on the data in a * single column.) Note that the output has P-values for * (a) the one-sample t-test, * (b) the sign test, and * (c) the Wilcoxon signed rank test * all for H_0:E(X)=0. In particular, we can carry out any of these * tests (or all of them) using proc univariate. * * Here we are mostly interested in the P-value for the Shapiro-Wilk * statistic W for normality. This test statistic has the unusual * property that SMALL VALUES suggest H_1 and thus rejection of * H_0, or departure from normality in this case. That is, the * P-value is from the LOWER tail: Pvalue=Prob(WWobs). * * The Shapiro-Wilk test statistic W is a ratio of two quantities. * The numerator has a `nonparametric' estimate of the variance, * which is divided by a sum-of-squares estimator of the variance. * * Small values of W suggest that the sums of squares overestimates * the variability of the data, so that the distribution of data * values is too `heavy-tailed' to be normally distributed. Since * the ANOVA test statistic depends heavily on sums of squares, this * tests for the most destructive way that data can deviate from * normality. * * A statistical test with a single P-value can, in principle, only * test for deviation as measured by a single one-dimensional * statistic. Since too-high a sum of squares is the most damaging, * the W test is the most appropriate test of normality for ANOVA * statistics. * * `Proc univariate' carries out tests for normality only if you say * `normal', as in `proc univariate normal'. Normal plots are * included if you add `plot'. These are plots with distorted * X and/or Y axes so that normally distributed data should appear * as a straight line. **********************************************************; proc univariate normal plot; title2 'PROC UNIVARIATE: Testing for normality'; title3 'The W (Shapiro-Wilk) test tests for normality'; title4 'The P-value is Prob(W(rand)