Power and Sample Size in SAS
Statistical Programs, University of Idaho
2022-05-03
1 Introduction
Estimation of sample size and statistical power rely on several components related to researcher’s objectives, the experimental design, the topic being studied, as well as the sampling methods used and the data being collected. Some of the statistical terms used in the process of estimation are described below.
n - Sample Size: This value is either estimated or specified by the researcher, as in the case of power estimation. It may be represented as either the size per treatment group, or the total sample size, depending on the estimation procedure used.
\(\sigma^2\) - Variance: This is a fixed value of variability for the data. The value used here may come from existing data, literature review, or simply the researcher’s best approximation. It is often expressed as a data standard deviation, \(s\), the square root of the mean squared error, \(\sqrt{MSE}\), or an estimated variance component from a mixed model.
\(\delta\) - Estimate of Precision: This is also a fixed value, set by the researcher, representing the degree of precision to be expected from the experimental estimates or the desired effect size to detect. It might be derived, for example, from the half width of a confidence interval, the expected or observed difference between treatment estimates, or a relative value such as a percentage of the estimate.
\(\alpha\) - Significance Level: The fixed level of significance for hypothesis testing, or conversely, (1- \(\alpha\)), the level of confidence coefficient desired. As defined by R. A. Fisher, \(\alpha\) is the probability of falsely rejecting the null hypothesis of no treatment effects (Type I error or false positive).
\(\beta\) - Type II Error Rate: Introduced into the statistical testing paradigm by statisticians Neyman and Pearson, \(\beta\) represents the probability of a false negative result. The converse probability, (1 – \(\beta\)), is the statistical power or the probability of correctly identifying a true rejection of the null hypothesis. Unlike the other components above, \(\beta\) is entirely dependent on the hypothesis being proposed. That is, the statistical power of any experiment will change depending on the hypothesis being tested. In sample size estimation, (1 - \(\beta\)) is specified by the researcher, while in power analysis, it is estimated from the other components and a fixed sample size.
In summary, \(\delta\), \(\alpha\), and \(\sigma\) are quantities predetermined by the researcher through known data, subject matter expertise, and expected outcomes. The sample size, n, and statistical power, (1 - \(\beta\)), are quantities to be estimated with the following restrictions: n can be estimated when \(\beta\) is specified, or conversely, (1 - \(\beta\)) can be estimated if n is given.
2 Basic Sample Size and Power Estimation
The SAS procedure Proc Power is designed to handle many different types of sample size and power computations. While this procedure can deal with more advanced designs such as ANOVA, linear regression and logistic regression, these techniques are now handled better through other procedures to be covered in following sections. In this tutorial, Proc Power will be used to address two types of problems: Estimation related to a single mean or proportion, and estimation/comparison of two independent means or proportions.
2.1 Single means
The Power procedure is called using standard SAS format with the Proc Power statement. In the example below for the case of a single mean, the statement ONESAMPLEMEANS is used with several options specifying the components outlined above. The first is the mean= option, giving the value anticipated for the mean estimate. The nullmean= option provides information relating to \(\delta\), the precision component. Several precision values can be specified at once in order to simultaneously generate several corresponding sample size or power estimates. In the example below, the option is called with the mean value set to 20, which will be evaluated with several precision values: 5, 10, 15 and 20. Note, however, the precisions are specified as the mean + \(\delta\) in the nullmean= option, e.g. 25, 30, 35, and 40.
The variability component, \(\sigma\), is given with the std= option at a value of 5 and the desired statistical power, 1 – \(\beta\), is fixed to 0.95. In this example, we are computing the estimated sample size, hence the option ntotal= is set to a missing value (a single period). This tells SAS that the sample size, n, is to be estimated. Alternatively, if we were computing estimated power, the ntotal= statement should be set to a predetermined value, while the power= option would be given a missing value, thereby informing SAS we want to estimate power. By default, \(\alpha\) is set to 0.05, but could be specified with an alpha= option.
proc power;
ONESAMPLEMEANS
mean=20
nullmean=25, 30, 35, 40
std=5
ntotal=.
power=.95;
run;
| Fixed Scenario Elements | |
|---|---|
| Distribution | Normal |
| Method | Exact |
| Mean | 20 |
| Standard Deviation | 5 |
| Nominal Power | 0.95 |
| Number of Sides | 2 |
| Alpha | 0.05 |
| Computed N Total | |||
|---|---|---|---|
| Index | Null Mean | Actual Power | N Total |
| 1 | 25 | 0.962 | 16 |
| 2 | 30 | 0.971 | 6 |
| 3 | 35 | 0.967 | 4 |
| 4 | 40 | 0.998 | 4 |
The output shows that, at the highest level of precision, \(\delta=5\), it is estimated to take at least 16 observations to reach 95% statistical power, while a larger precision of \(\delta=20\) would require only 4 observations. The Actual Power is an estimate of (1 - \(\beta\)) for each sample size. As a general rule, the smaller the desired precision is, the larger the sample size will need to be.
2.2 Single proportions
Tests for proportions can be run in a similar manner using the ONESAMPLEFREQ statement and corresponding proportion = and nullproportion= options to specify the expected proportion and precision(s), respectively.
In the example below, we now want to estimate power given a sample size. It is possible here to specify several sample sizes. In addition, we also specify several precision values in the nullproportion= option using range notation, i.e. <lower value> to <upper value> by <increment>. This notation is a shorthand syntax so we don’t have to type out each value separately. The expected proportion in this example is set to 0.31 or 12/39. The sides=1 option specifies a one-sided hypothesis test. That is, we are only interested in values greater than 0.31.
Because we are requesting multiple power estimates with all these settings, it is often easier to review the results graphically in a Power Curve. The ODS statement and code following Proc Power will set up this plot.
proc power;
onesamplefreq
proportion=0.31
nullproportion = 0.01 to 0.3 by 0.01
ntotal = 20, 30, 50
sides=1
power=.;
ods output Output=power;
run;
proc sort data=power;
by nullproportion NTotal ;
run;
proc transpose data=power out=power;
var power;
id NTotal;
by nullproportion;
run;
data power;
set power;
delta = 12/39 - nullproportion;
run;
proc sgplot;
series x=delta y = _20/name='twenty' legendlabel="N=20";
series x=delta y = _30/ legendlabel="N=30";
series x=delta y = _50/ legendlabel="N=50";
refline 0.8 /axis=Y lineattrs=(pattern=shortdash color=black);
xaxis label='Detectable difference at 95% confidence' LABELATTRS=(Family=Arial Size=1 Weight=Bold) VALUEATTRS=(Family=Arial Size=12 Weight=Bold);
yaxis label='Power' min=0 LABELATTRS=( Family=Arial Size=15 Weight=Bold) VALUEATTRS=(Family=Arial Size=12 Weight=Bold);
run;
| Fixed Scenario Elements | |
|---|---|
| Method | Exact |
| Number of Sides | 1 |
| Binomial Proportion | 0.31 |
| Nominal Alpha | 0.05 |
| Computed Power | ||||||
|---|---|---|---|---|---|---|
| Index | Null Proportion | N Total |
Lower Crit Val |
Upper Crit Val |
Actual Alpha | Power |
| 1 | 0.01 | 20 | . | 2 | 0.01686 | 0.994 |
| 2 | 0.01 | 30 | . | 2 | 0.03615 | >.999 |
| 3 | 0.01 | 50 | . | 3 | 0.01382 | >.999 |
| 4 | 0.02 | 20 | . | 3 | 0.00707 | 0.971 |
| 5 | 0.02 | 30 | . | 3 | 0.02172 | 0.999 |
| 6 | 0.02 | 50 | . | 4 | 0.01776 | >.999 |
| 7 | 0.03 | 20 | . | 3 | 0.02101 | 0.971 |
| 8 | 0.03 | 30 | . | 4 | 0.01190 | 0.993 |
| 9 | 0.03 | 50 | . | 5 | 0.01681 | >.999 |
| 10 | 0.04 | 20 | . | 3 | 0.04386 | 0.971 |
| 11 | 0.04 | 30 | . | 4 | 0.03059 | 0.993 |
| 12 | 0.04 | 50 | . | 5 | 0.04897 | >.999 |
| 13 | 0.05 | 20 | . | 4 | 0.01590 | 0.909 |
| 14 | 0.05 | 30 | . | 5 | 0.01564 | 0.977 |
| 15 | 0.05 | 50 | . | 6 | 0.03778 | >.999 |
| 16 | 0.06 | 20 | . | 4 | 0.02897 | 0.909 |
| 17 | 0.06 | 30 | . | 5 | 0.03154 | 0.977 |
| 18 | 0.06 | 50 | . | 7 | 0.02892 | 0.998 |
| 19 | 0.07 | 20 | . | 4 | 0.04713 | 0.909 |
| 20 | 0.07 | 30 | . | 6 | 0.01623 | 0.939 |
| 21 | 0.07 | 50 | . | 8 | 0.02201 | 0.995 |
| 22 | 0.08 | 20 | . | 5 | 0.01834 | 0.791 |
| 23 | 0.08 | 30 | . | 6 | 0.02929 | 0.939 |
| 24 | 0.08 | 50 | . | 8 | 0.04379 | 0.995 |
| 25 | 0.09 | 20 | . | 5 | 0.02904 | 0.791 |
| 26 | 0.09 | 30 | . | 6 | 0.04806 | 0.939 |
| 27 | 0.09 | 50 | . | 9 | 0.03283 | 0.987 |
| 28 | 0.10 | 20 | . | 5 | 0.04317 | 0.791 |
| 29 | 0.10 | 30 | . | 7 | 0.02583 | 0.867 |
| 30 | 0.10 | 50 | . | 10 | 0.02454 | 0.971 |
| 31 | 0.11 | 20 | . | 6 | 0.01755 | 0.621 |
| 32 | 0.11 | 30 | . | 7 | 0.04074 | 0.867 |
| 33 | 0.11 | 50 | . | 10 | 0.04345 | 0.971 |
| 34 | 0.12 | 20 | . | 6 | 0.02602 | 0.621 |
| 35 | 0.12 | 30 | . | 8 | 0.02205 | 0.757 |
| 36 | 0.12 | 50 | . | 11 | 0.03250 | 0.941 |
| 37 | 0.13 | 20 | . | 6 | 0.03697 | 0.621 |
| 38 | 0.13 | 30 | . | 8 | 0.03390 | 0.757 |
| 39 | 0.13 | 50 | . | 12 | 0.02419 | 0.892 |
| 40 | 0.14 | 20 | . | 7 | 0.01534 | 0.431 |
| 41 | 0.14 | 30 | . | 8 | 0.04967 | 0.757 |
| 42 | 0.14 | 50 | . | 12 | 0.04022 | 0.892 |
| 43 | 0.15 | 20 | . | 7 | 0.02194 | 0.431 |
| 44 | 0.15 | 30 | . | 9 | 0.02778 | 0.615 |
| 45 | 0.15 | 50 | . | 13 | 0.03006 | 0.820 |
| 46 | 0.16 | 20 | . | 7 | 0.03037 | 0.431 |
| 47 | 0.16 | 30 | . | 9 | 0.04028 | 0.615 |
| 48 | 0.16 | 50 | . | 13 | 0.04755 | 0.820 |
| 49 | 0.17 | 20 | . | 7 | 0.04088 | 0.431 |
| 50 | 0.17 | 30 | . | 10 | 0.02243 | 0.459 |
| 51 | 0.17 | 50 | . | 14 | 0.03573 | 0.725 |
| 52 | 0.18 | 20 | . | 8 | 0.01771 | 0.259 |
| 53 | 0.18 | 30 | . | 10 | 0.03231 | 0.459 |
| 54 | 0.18 | 50 | . | 15 | 0.02664 | 0.613 |
| 55 | 0.19 | 20 | . | 8 | 0.02413 | 0.259 |
| 56 | 0.19 | 30 | . | 10 | 0.04508 | 0.459 |
| 57 | 0.19 | 50 | . | 15 | 0.04111 | 0.613 |
| 58 | 0.20 | 20 | . | 8 | 0.03214 | 0.259 |
| 59 | 0.20 | 30 | . | 11 | 0.02562 | 0.311 |
| 60 | 0.20 | 50 | . | 16 | 0.03080 | 0.492 |
| 61 | 0.21 | 20 | . | 8 | 0.04192 | 0.259 |
| 62 | 0.21 | 30 | . | 11 | 0.03571 | 0.311 |
| 63 | 0.21 | 50 | . | 16 | 0.04615 | 0.492 |
| 64 | 0.22 | 20 | . | 9 | 0.01859 | 0.134 |
| 65 | 0.22 | 30 | . | 11 | 0.04848 | 0.311 |
| 66 | 0.22 | 50 | . | 17 | 0.03474 | 0.373 |
| 67 | 0.23 | 20 | . | 9 | 0.02461 | 0.134 |
| 68 | 0.23 | 30 | . | 12 | 0.02795 | 0.191 |
| 69 | 0.23 | 50 | . | 18 | 0.02589 | 0.267 |
| 70 | 0.24 | 20 | . | 9 | 0.03200 | 0.134 |
| 71 | 0.24 | 30 | . | 12 | 0.03806 | 0.191 |
| 72 | 0.24 | 50 | . | 18 | 0.03843 | 0.267 |
| 73 | 0.25 | 20 | . | 9 | 0.04093 | 0.134 |
| 74 | 0.25 | 30 | . | 13 | 0.02159 | 0.105 |
| 75 | 0.25 | 50 | . | 19 | 0.02873 | 0.179 |
| 76 | 0.26 | 20 | . | 10 | 0.01831 | 0.059 |
| 77 | 0.26 | 30 | . | 13 | 0.02950 | 0.105 |
| 78 | 0.26 | 50 | . | 19 | 0.04183 | 0.179 |
| 79 | 0.27 | 20 | . | 10 | 0.02379 | 0.059 |
| 80 | 0.27 | 30 | . | 13 | 0.03947 | 0.105 |
| 81 | 0.27 | 50 | . | 20 | 0.03138 | 0.112 |
| 82 | 0.28 | 20 | . | 10 | 0.03047 | 0.059 |
| 83 | 0.28 | 30 | . | 14 | 0.02254 | 0.052 |
| 84 | 0.28 | 50 | . | 20 | 0.04495 | 0.112 |
| 85 | 0.29 | 20 | . | 10 | 0.03847 | 0.059 |
| 86 | 0.29 | 30 | . | 14 | 0.03033 | 0.052 |
| 87 | 0.29 | 50 | . | 21 | 0.03380 | 0.066 |
| 88 | 0.30 | 20 | . | 10 | 0.04796 | 0.059 |
| 89 | 0.30 | 30 | . | 14 | 0.04005 | 0.052 |
| 90 | 0.30 | 50 | . | 21 | 0.04776 | 0.066 |
As can be seen, the tabular output from all these options can be verbose. The Power Curve plot, however, summarizes things neatly. We can see that with a sample size of n=50, we reach 80% power with a precision of about \(\delta=0.15\), while similar power for n=20 will only give a precision of about \(\delta=0.25\).
2.3 Comparison of Two means or proportions
As in the single estimate case, there are different statements for comparing two means or proportions: TWOSAMPLEMEANS and TWOSAMPLEFREQ, respectively.
Options in these statements now provide information on pairs or groups. For the example below, the groupmeans= option specifies the means for two samples. They are separated by the vertical line character ” | “. Here, one mean is set to 20, while the second mean is evaluated at four values: 25, 30, 35, and 40. While this is similar to the one sample case above, we are now estimating the sample size for the difference of two means, rather than the sample size for estimating one mean at varying precision values. That is, the sample size will be evaluated for four tests: 20 vs 25, 20 vs 30, 20 vs 35, and 20 vs 40. The precision component, \(\delta\), in the two sample problem is now represented as the difference between the two estimated quantities (effect size). Similar to the means, the variability must also now be specified for two groups using the groupstds= option. In this case, they are set to 5 and 7. The test=diff_satt tells SAS that the hypothesis to be tested is the difference of the two means and that the standard deviations of each group could be different. The sample size can be specified using either the ntotal= or npergroup= option. The power= option is specified as before.
proc power;
twosamplemeans
groupmeans=20|25 30 35 40
groupstds=5|7
Test=DIFF_SATT
npergroup=.
power=.95;
run;
| Fixed Scenario Elements | |
|---|---|
| Distribution | Normal |
| Method | Exact |
| Group 1 Mean | 20 |
| Group 1 Standard Deviation | 5 |
| Group 2 Standard Deviation | 7 |
| Nominal Power | 0.95 |
| Number of Sides | 2 |
| Null Difference | 0 |
| Nominal Alpha | 0.05 |
| Computed N per Group | ||||
|---|---|---|---|---|
| Index | Mean2 | Actual Alpha | Actual Power | N per Group |
| 1 | 25 | 0.0500 | 0.952 | 40 |
| 2 | 30 | 0.0493 | 0.953 | 11 |
| 3 | 35 | 0.0471 | 0.961 | 6 |
| 4 | 40 | 0.0457 | 0.990 | 5 |
The effect of precision can easily be seen here. For a high precision of \(\delta = 5\), or \(\mu_1 = 20\) and \(\mu_2 = 25\), we need 40 observations per group, while a precision of \(\delta = 20\) requires only 5 observations per group.
3 Advanced Sample Size and Power Estimation
In this section, power and sample size issues are addressed for more sophisticated experimental designs and analyses. The examples below will deal with situations involving analysis of variance, although the similar examples can typically be adapted to cases for linear and multiple regressions as well.
3.1 Evaluating Existing Data
At times, we wish to examine the statistical power for an experiment that has already been carried out. In this first example, the randomized complete block (RCB) experimental design is considered using data previously shown in the Proc Mixed tutorial here.
These data describe a variety (2 levels) and fertilizer (3 levels) study with 4 blocks measuring yield. We will take the estimated means from this analysis:
| Variety | Fertilizer | Yield LSMean | |
|---|---|---|---|
| VAR1 | 1 | 1.85871031 | |
| VAR1 | 2 | 3.12521718 | |
| VAR1 | 3 | 3.92756253 | |
| VAR2 | 1 | 0.86258303 | |
| VAR2 | 2 | 2.94528678 | |
| VAR2 | 3 | 3.17151761 |
as well as the estimated variance, \(s^2 = 0.6077\) to compute power for this data under various sample size scenarios. Because of the factorial treatment design (2 varieties x 3 fertilizers), there will be multiple hypotheses to compute power for: Variety, Fertilizer, and the interaction Variety*Fertilizer. In addition, we will add contrasts that compare specific levels of the significant effect, Fertilizer. The SAS procedure for these power calculations is Proc GLMPower. This procedure is based on a fixed linear model, not a mixed linear model, however, it can still provide some information on sample sizes and power. Specific procedures for mixed model estimation are given in subsequent sections.
We begin the analysis by entering the estimated means into a data set. The GLMPower procedure is then called.
The Class, Model and Contrast statements are similar to those that would be used in Proc Mixed. The Power statement specifies necessary information for power calculations. For example, the stddev is given as the square root of the residual variance: \(\sqrt{0.6077} = 0.7795\). There were originally 2x3 = 6 treatment combinations at 4 blocks each giving a total sample size of ntotal = 24. The power is set to missing in order to solve for that component. The probability of Type I error rate, \(\alpha\), is assumed to be 0.05 and the precision, \(\delta\), is determined from the means in the data set. The plot statement produces a power plot of all the specified hypotheses. The statistical power is displayed on the Y axis and sample size, \(n\), on the X axis. The X axis has been specified to vary from \(n = 24\) to \(n = 275\).
data yield_means;
input var$ fert yield;
cards;
VAR1 1 1.8587103
VAR1 2 3.1252172
VAR1 3 3.9275625
VAR2 1 0.8625830
VAR2 2 2.9452868
VAR2 3 .1715176
;
run;
proc glmpower data=yield_means;
class var fert;
model yield = var fert var*fert;
contrast 'Fert 1 vs 3' fert 1 0 -1;
contrast 'Fert 2 vs 3' fert 0 1 -1;
power
stddev = 0.7795
ntotal = 24
power = .;
plot x=n min=24 max=275;
title1 'GLMPower Example';
run;| GLMPower Example |
| Fixed Scenario Elements | |
|---|---|
| Dependent Variable | yield |
| Error Standard Deviation | 0.7795 |
| Total Sample Size | 24 |
| Alpha | 0.05 |
| Error Degrees of Freedom | 18 |
| Computed Power | ||||
|---|---|---|---|---|
| Index | Type | Source | Test DF | Power |
| 1 | Effect | var | 1 | 0.998 |
| 2 | Effect | fert | 2 | 0.953 |
| 3 | Effect | var*fert | 2 | 0.982 |
| 4 | Contrast | Fert 1 vs 3 | 1 | 0.387 |
| 5 | Contrast | Fert 2 vs 3 | 1 | 0.667 |
| GLMPower Example |
The table indicates that there are limited powers for the Variety and Variety*Fertilizer interaction hypotheses (\(0.482\) and \(0.131\), respectively). The Fertilizer effect, however, has very high power \((0.999)\). Comparison of Fertilizer means shows that there is high power for Fertilizer 1 vs Fertilizer 3 \((0.999)\) and limited power for Fertilizer 2 vs Fertilizer 3 \((0.239)\).
Examining these hypotheses across sample sizes in the corresponding plot shows flat lines for the overall Fertilizer effect and the contrast of Fertilizer 1 vs Fertilizer 3 (top of the figure). These hypotheses have very high power (ability to detect true differences) at all sample sizes. The Variety effect (solid line) has low power at \(n = 24 (0.482)\) and does not obtain reasonable power until \(n ≥ 50\), that is, at least 8 blocks. The last two hypotheses, the Variety*Fertilizer interaction and the Fertilizer 2 vs Fertilizer 3 contrast, do not have good power with less than \(n = 150\) or 26 blocks. Hence, we can conclude that this experiment only had good power for large differences in the fertilizer effect and that very large sample sizes would have been required in order to detect other effects.
3.2 Mixed Models
In most experimental settings, it is often more appropriate to specify a mixed model possibly with a non-normal response in order to fully account for the underlying data structure. For example, while the normality and independence assumptions used in the RCB demonstration shown above were appropriate, the effect of blocks would be better represented as a random effect. In such cases, another means of assessing statistical power is required. The RCB example can be reanalyzed in a mixed model format using the following PROC GLIMMIX code:
proc glimmix data=yield1;
class block var fert;
model yield = var fert var*fert/dist=normal;
random block;
run;where, block is now assumed to be a random effect (for more information on generalized linear mixed models, please refer to the GLMMIX tutorial.
In this analysis, the estimates for the means and residual variance are the same as given above. We also have a variance component estimate for block of 0.5457. This information can then be used to estimate statistical power from a technique outlined by Stroup (1999) that utilizes the theoretical relationship between the F statistic, degrees of freedom, and a value called the non-centrality parameter (NC). The steps to do this are given below.
The analysis means are put into a dataset as before with the addition of Do loop to create a block variable with values 1, 2, 3, and 4. This essentially just repeats the means once for each block 4 times.
Proc Glimmix is called as above with the addition of an option noprofile. Two additional statements are also present. The parms statement lists the values for the random effects (Blocks and the residual error, in that order). The hold=1,2 option requests that parameters 1 and 2 be held constant, telling SAS to just use these values and not to estimate them for the data. The ODS statement outputs the test statistics to a data set “results” for further computations. Those computations are listed in the data step that followsand involves the theoretical relationship between the F statistic, degrees of freedom, and the non-centrality parameter (NC) mentioned earlier. For the user, all that needs to be specified in this data step is the value for alpha. Lastly, the computed power data set is printed.
data rcb_means;
input var$ fert yield;
do block = 1 to 4;
output;
end;
cards;
VAR1 1 1.8587103
VAR1 2 3.1252172
VAR1 3 3.9275625
VAR2 1 0.8625830
VAR2 2 2.9452868
VAR2 3 3.1715176
run;
proc glimmix data=rcb_means noprofile;
class block var fert;
model yield= var fert var*fert/dist=normal;
random block;
parms (0.5457) (0.6077)/hold=1,2;
ods output tests3=results;
run;
data power;
set results;
alpha = 0.05;
nc = numdf*fvalue;
fcrit = finv(1 - alpha, numdf, dendf, 0);
power = 1 - probf(fcrit, numdf, dendf, nc);
run;
proc print;
run;
| Model Information | |
|---|---|
| Data Set | WORK.RCB_MEANS |
| Response Variable | yield |
| Response Distribution | Gaussian |
| Link Function | Identity |
| Variance Function | Default |
| Variance Matrix | Not blocked |
| Estimation Technique | Restricted Maximum Likelihood |
| Degrees of Freedom Method | Containment |
| Class Level Information | ||
|---|---|---|
| Class | Levels | Values |
| block | 4 | 1 2 3 4 |
| var | 2 | VAR1 VAR2 |
| fert | 3 | 1 2 3 |
| Number of Observations Read | 24 |
|---|---|
| Number of Observations Used | 24 |
| Dimensions | |
|---|---|
| G-side Cov. Parameters | 1 |
| R-side Cov. Parameters | 1 |
| Columns in X | 12 |
| Columns in Z | 4 |
| Subjects (Blocks in V) | 1 |
| Max Obs per Subject | 24 |
| Parameter Search | ||
|---|---|---|
| CovP1 | CovP2 |
Objective Function |
| 0.5457 | 0.6077 | 37.997418438 |
| Optimization Information | |
|---|---|
| Optimization Technique | Dual Quasi-Newton |
| Parameters in Optimization | 2 |
| Equality Constraints | 2 |
| Lower Boundaries | 2 |
| Upper Boundaries | 2 |
| Fixed Effects | Profiled |
| Starting From | Data |
| Iteration History | |||||
|---|---|---|---|---|---|
| Iteration | Restarts | Evaluations |
Objective Function |
Change |
Max Gradient |
| 0 | 0 | 4 | 37.997418438 | . | 0 |
| Convergence criterion (ABSGCONV=0.00001) satisfied. |
| Fit Statistics | |
|---|---|
| -2 Res Log Likelihood | 38.00 |
| AIC (smaller is better) | 38.00 |
| AICC (smaller is better) | 38.00 |
| BIC (smaller is better) | 38.00 |
| CAIC (smaller is better) | 38.00 |
| HQIC (smaller is better) | 38.00 |
| Generalized Chi-Square | -0.00 |
| Gener. Chi-Square / DF | -0.00 |
| Covariance Parameter Estimates | ||
|---|---|---|
| Cov Parm | Estimate |
Standard Error |
| block | 0.5457 | . |
| Residual | 0.6077 | . |
| Type III Tests of Fixed Effects | ||||
|---|---|---|---|---|
| Effect | Num DF | Den DF | F Value | Pr > F |
| var | 1 | 15 | 4.10 | 0.0612 |
| fert | 2 | 15 | 17.25 | 0.0001 |
| var*fert | 2 | 15 | 0.58 | 0.5725 |
| Obs | Effect | NumDF | DenDF | FValue | ProbF | alpha | nc | fcrit | power |
|---|---|---|---|---|---|---|---|---|---|
| 1 | var | 1 | 15 | 4.10 | 0.0612 | 0.05 | 4.0952 | 4.54308 | 0.47365 |
| 2 | fert | 2 | 15 | 17.25 | 0.0001 | 0.05 | 34.4909 | 3.68232 | 0.99844 |
| 3 | var*fert | 2 | 15 | 0.58 | 0.5725 | 0.05 | 1.1582 | 3.68232 | 0.12847 |
The power estimates are very similar to those seen previously, but differ slightly due to differences in the estimation procedure and the addition of a random effect for blocks.
3.3 Mixed Models with Non-normal Responses
The real value of a procedure such as Proc Glimmix, however, is its ability to handle non-normal responses. The following example is, again, a RCB design with two treatments (0 and 1) and 8 blocks. The response in this case, however, is binomial, measuring a positive or negative result. Each treatment was recorded 100 times per block. An analysis of variance assuming a generalized linear mixed model for a binomial response and a logit link function was carried out for these data and can be seen in the Proc Glimmix tutorial here.
As was shown above, we start by using the estimated means (inverse link values) and the random effects estimates for block and block*treatment. The means are put into a dataset and this time, the number of blocks is allowed to range from 2 to 16. Within each block we fix the total number of binomial trials, N, in each treatment and block to 200 and compute an expected number of counts, Y, as this total x the estimated mean proportion. The subsequent power computation is then carried out and results are printed and plotted.
data intro_binomial;
input Treatment P;
N = 200;
Y = N*P;
do B = 2 to 16;
do block = 1 to B;
output;
end;
end;
cards;
0 0.9276
1 0.7807
;
run;
proc sort;
by B;
proc glimmix data=intro_binomial noprofile;
class block Treatment;
model Y/N= Treatment;
random block block*Treatment;
parms (0.5201) (0.8335)/hold=1,2;
ods output tests3=results;
by B;
run;
data power;
set results;
alpha = 0.05;
nc = numdf*fvalue;
fcrit = finv(1 - alpha, numdf, dendf, 0);
power = 1 - probf(fcrit, numdf, dendf, nc);
run;
proc print;
run;
proc sgplot;
scatter x=b y=power;
series x=b y=power;
refline .80/axis=y;
xaxis label='Number of blocks' LABELATTRS=(Family=Arial Size=13 Weight=Bold) VALUEATTRS=(Family=Arial Size=12 Weight=Bold);
yaxis label='Power' min=0 LABELATTRS=(Family=Arial Size=15 Weight=Bold) VALUEATTRS=(Family=Arial Size=12 Weight=Bold);
run;
| Obs | B | Effect | NumDF | DenDF | FValue | ProbF | alpha | nc | fcrit | power |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2 | Treatment | 1 | 1 | 1.85 | 0.4034 | 0.05 | 1.8525 | 161.448 | 0.09005 |
| 2 | 3 | Treatment | 1 | 2 | 2.78 | 0.2375 | 0.05 | 2.7787 | 18.513 | 0.17035 |
| 3 | 4 | Treatment | 1 | 3 | 3.70 | 0.1499 | 0.05 | 3.7050 | 10.128 | 0.27210 |
| 4 | 5 | Treatment | 1 | 4 | 4.63 | 0.0978 | 0.05 | 4.6312 | 7.709 | 0.37772 |
| 5 | 6 | Treatment | 1 | 5 | 5.56 | 0.0650 | 0.05 | 5.5574 | 6.608 | 0.47804 |
| 6 | 7 | Treatment | 1 | 6 | 6.48 | 0.0437 | 0.05 | 6.4837 | 5.987 | 0.56874 |
| 7 | 8 | Treatment | 1 | 7 | 7.41 | 0.0297 | 0.05 | 7.4099 | 5.591 | 0.64811 |
| 8 | 9 | Treatment | 1 | 8 | 8.34 | 0.0203 | 0.05 | 8.3361 | 5.318 | 0.71594 |
| 9 | 10 | Treatment | 1 | 9 | 9.26 | 0.0139 | 0.05 | 9.2624 | 5.117 | 0.77284 |
| 10 | 11 | Treatment | 1 | 10 | 10.19 | 0.0096 | 0.05 | 10.1886 | 4.965 | 0.81984 |
| 11 | 12 | Treatment | 1 | 11 | 11.11 | 0.0067 | 0.05 | 11.1149 | 4.844 | 0.85818 |
| 12 | 15 | Treatment | 1 | 14 | 13.89 | 0.0023 | 0.05 | 13.8936 | 4.600 | 0.93341 |
| 13 | 16 | Treatment | 1 | 15 | 14.82 | 0.0016 | 0.05 | 14.8198 | 4.543 | 0.94881 |
This shows that the number of blocks needs to be 11 or more to reach a statistical power for the treatment effect of at least 0.80.
This technique can be useful when dealing with non-normal data, or situations where there are complex model and covariance structures. The user will, however, need to have some estimates or likely values for all means and variance components.