SEPTEMBER2007
DOI: 10.1007/S11336-007-9012-5
A UNIFIED APPROACH TO POWER CALCULATION AND SAMPLE SIZE DETERMINATION FOR RANDOM REGRESSION MODELS
G
WOWENS
HIEHNATIONAL CHIAO TUNG UNIVERSITY
The underlying statistical models for multiple regression analysis are typically attributed to two types of modeling: fixed and random. The procedures for calculating power and sample size under the fixed regression models are well known. However, the literature on random regression models is limited and has been confined to the case of all variables having a joint multivariate normal distribution. This paper presents a unified approach to determining power and sample size for random regression models with arbitrary distribution configurations for explanatory variables. Numerical examples are provided to illustrate the usefulness of the proposed method and Monte Carlo simulation studies are also conducted to assess the accuracy. The results show that the proposed method performs well for various model speci-fications and explanatory variable distributions.
Key words: asymptotic distribution, effect size, noncentral F distribution.
1. Introduction
Multiple regression analysis is one of the widely used statistical methods. Conventionally, there are two approaches to the statistical modeling of these regression applications. They are referred to as fixed (conditional) and random (unconditional) models. In the context of regres-sion analysis, it is quite common in the behavioral and social sciences to have studies in which not only the values of response variables for each experimental unit are just available after the observations are made, but also the levels of explanatory variables cannot be fixed in advance. Therefore, the explanatory variables are also outcomes of the study under such circumstances. In order to take account of this extra variability, the appropriate strategy is to consider the random regression setting. On the other hand, the fixed regression model is suitable for studies in which the configurations of the explanatory variables are preset by the researcher.
Sample size calculations and power analyses are often critical for researchers to address spe-cific scientific hypotheses and confirm credible treatment effects. Thus, they should be an integral part of the whole study. Accordingly, it is of practical importance to be able to perform these tasks in a multiple regression setup. For fixed regression models, the procedures are well documented in the literature. Regarding random regression models, Gatsonis and Sampson (1989) gave an ex-cellent and thorough description of exact power and sample size calculations when the response and explanatory variables have a joint multivariate normal distribution. Traditionally, the problem is referred to as multiple correlation analysis and the parameter of interest is the squared multiple correlation coefficient. In contrast, Algina and Olejnik (2000) presented results for determining the sample size required for adequate estimation accuracy of the squared multiple correlation co-efficient. Furthermore, related treatments can be found in Kelley and Maxwell (2003), Mendoza and Stafford (2001), Shieh (2006) and Steiger and Fouladi (1992). It is important to note that
The author would like to thank the editor, the associate editor, and the referees for drawing attention to pertinent references that led to improved presentation. This research was partially supported by National Science Council grant NSC-94-2118-M-009-004.
Requests for reprints should be sent to Gwowen Shieh, Department of Management Science, National Chiao Tung University, Hsinchu, Taiwan 30050, ROC. E-mail: gwshieh@mail.nctu.edu.tw
the studies cited above for random regression models and multiple correlation analysis are ap-plicable in the circumstance that the response and explanatory variables have a joint multivariate normal distribution. However, there are many situations in which assuming normal distribution for explanatory variables is inappropriate. For instance, consider the simple interaction model in the formulation of Y = β0+ Xβ1+ Zβ2+ XZβ3+ ε. It is commonly assumed that the con-tinuous measurements X and Z are normally distributed. However, the product of two normally distributed variables (XZ) does not have a normal distribution. Therefore, the existing results for power analysis of a multinormal situation do not apply in this application. In fact, Gatsonis and Sampson (1989) noted that when the joint distribution of Y and all explanatory variables is non-normal, it is doubtful that their power calculations give results that are accurate and reasonable. Therefore, neglecting other configurations of explanatory variables is an obvious limitation of available methods. A natural generalization to incorporate both normal and nonnormal explana-tory variables should be essential to the existing approaches for performing power and sample size calculations in practice. It should be observed that the prescribed simple interaction model is directly connected to the moderated multiple regression formulation which has been pervasively used for testing moderator effects in all areas of social sciences, see Aguinis, Beaty, Boik, and Pierce (2005) for further details.
This paper aims to provide a unified approach to the determinations of power and sample size for random regression models. The distinct feature of the proposed method is the accommo-dation of arbitrary discrete and/or continuous distribution formulations for explanatory variables. Therefore, the aforementioned multivariate normal setting can be viewed as a special case. For related results and various extensions in hierarchical linear models and multivariate linear mod-els, the interested reader is referred to Raudenbush and Liu (2000,2001), Shieh (2003,2005), and the references therein. In fact, the suggested two-stage methodology can be viewed as an extension of the results for the univariate case in Shieh (2005). The rest of the paper is organized as follows. In Section2, the important analytical details of the proposed method are described. Numerical examples are provided in Section3to demonstrate the proposed power and sample size calculations for several random regression formulations. Since the approach considered here uses large sample approximations, simulation studies are conducted to assess its adequacy for finite sample and robustness under various model specifications and distributions of explanatory variables. Finally, Section4contains some final remarks.
2. The Proposed Method
To facilitate the illustration of the proposed method for random regression models, it is instructive to review first the situation under the fixed regression models where the results would be specific to the particular values of the explanatory variables that are observed or predetermined by the researcher.
2.1. Review of Fixed Regression Models
Consider the standard multiple linear regression model with response variable Y and all the levels of p explanatory variables X(1), . . . , X(p) fixed a priori:
Y= Xβ + ε, (1)
where Y= (Y1, . . . , YN)T, Yi is the value of the response variable Y ; X= (1N,XD) where 1N is the N × 1 vector of all 1’s, XD = (X1, . . . ,XN)T is often called the design matrix, Xi = (xi1, . . . , xip)T, xi1, . . . , xip are the known constants of the p explanatory variables for i= 1, . . . , N; β = (β0, β1, . . . , βp)T where β0, β1, . . . , βp are unknown parameters; and ε=
(ε1, . . . , εN)Twhere εi are iid N (0, σ2)random variables. We are concerned with the general linear hypothesis H0: Lβ= θ versus H1: Lβ= θ, where L is an l × (p + 1) coefficient matrix of rank l≤ p + 1 and θ is an l × 1 vector of constants. It is well known that under the assumption given in (1), the likelihood ratio test for H0is based on
F = SSH/ l
SSE/(N− p − 1), (2)
where SSH= (L ˆβ − θ)T[L(XTX)−1LT]−1(L ˆβ− θ), SSE = (Y − X ˆβ)T(Y− X ˆβ) and ˆβ = (XTX)−1XTY is an unbiased estimator of β, see Rencher (2000, Chaps. 7–8) for further de-tails. Under the alternative hypothesis, F is distributed as F (l, N− p − 1, λ), the noncentral F distribution with l and N− p − 1 degrees of freedom and noncentrality parameter
λ= (Lβ − θ)TL(XTX)−1LT−1(Lβ− θ)/σ2. (3) If the null hypothesis is true, then λ= 0 and F is distributed as F (l, N − p − 1), a central or regular F distribution with l and N− p − 1 degrees of freedom. The test is carried out by rejecting H0if F > Fl,N−p−1,α, where Fl,N−p−1,α is the upper 100α percentage point of the central F distribution F (l, N− p − 1).
To calculate power and sample size, it is assumed that there are m distinct configurations of Xi for i= 1, . . . , N, and they are denoted by Zj with the proportions wj, j= 1, . . . , m (≤ N). Then, XTX can be expressed as XTX= N · , where =m
j=1wjZjZTj. Accordingly, the noncentrality parameter λ in (3) is rewritten as
λ= Nδ, (4)
where δ= (Lβ − θ)T[L−1LT]−1(Lβ− θ)/σ2is the so-called effect size, see Cohen (1988). Hence, given all model configurations and sample size N , the statistical power achieved for test-ing hypothesis H0: Lβ= θ with specified significance level α against the alternative H1: Lβ= θ is the probability
PF (l, N− p − 1, Nδ) > Fl,N−p−1,α
, (5)
where δ is defined in (4). Furthermore, this power function can be utilized to calculate the sample size needed in order to attain the specified power. However, it usually involves an iterative process to find the solution because both F (l, N− p − 1, Nδ) and Fl,N−p−1,α depend on the sample size N .
2.2. Random Regression Models
To extend the concept and interpretation of the aforementioned results, we assume that the explanatory variables{X∗i = x∗i, i= 1, . . . , N} in (1) have a probability function f (X∗i)with finite moments. The form of f (X∗i)is assumed to be dependent on none of the unknown para-meters β and σ2. Thus, the notations of Xi and XD as observed values in the fixed regression model are replaced by X∗i and X∗D= (X∗1, . . . ,X∗N)Tas random variables hereinafter. Frequently, the inferences are concerned mainly with the regression coefficients β1= (β1, . . . , βp)Tand the corresponding coefficient matrix is written in the form of L= L1, where L1= (0c,C), 0c is the c× 1 null vector of all 0’s and C is a c × p coefficient matrix of rank c ≤ p. It follows from the overall estimator ˆβ given above that the prescribed estimator for β1can be expressed as ˆβ∗1= (X∗TCX∗C)−1XC∗TY, where X∗C= (IN− J/N)X∗D is the centered form of X∗D, IN is the
identity matrix of dimension N and J is the N× N square matrix of all 1’s. In view of the extra random nature of X∗D, it is easily seen that
C ˆβ∗1|X∗D∼ Nc
Cβ1, σ2C
X∗TC X∗C−1CT.
It therefore follows that the general linear hypothesis reduces to H0: Cβ1= θ versus H1: Cβ1= θand the test statistic is of the form
F∗= SSH
∗/c
SSE∗/(N− p − 1), (6)
SSH∗= (C ˆβ∗1− θ)T[C(XC∗TX∗C)−1CT]−1(C ˆβ∗1− θ), SSE∗= (Y − X∗ˆβ∗)T(Y− X∗ˆβ∗), X∗= (1N,X∗D)and ˆβ∗= (X∗TX∗)−1X∗TY. Under random formulation, F∗has the conditional distri-bution
F∗| ∼ F (c, N − p − 1, ), (7)
where the noncentrality parameter = (Cβ1− θ)T[C(X∗TCX∗C)−1CT]−1(Cβ1− θ)/σ2is a ran-dom variable and the exact distribution depends ultimately on the joint distribution of X∗D. In order to provide a generally useful and versatile solution without specifically confining to any particular X∗D, the asymptotic property of is studied next.
Let μ and denote the mean vector and covariance matrix of the random explanatory vari-ables X∗i = (Xi1∗, . . . , X∗ip)T, respectively. It follows from the standard asymptotic result (Muir-head,1982, Cor. 1.2.18) that S∗= (X∗TCX∗C)/(N− 1) has asymptotic normal distribution
(N− 1)1/2vec(S∗)− vec() .∼ Np2
0p2, − vec() · vec()T
,
where vec(·) is a matrix operator which arranges the columns of a matrix into one long column,
= E[(X∗i − μ)(X∗i − μ)T⊗ (X∗i − μ)(X∗i − μ)T], E[·] denotes the expectation taken with
respect to the distribution of X∗i, and⊗ represents the Kronecker product. Using the identity vec(ABC)= (CT⊗ A) · vec(B), the noncentrality parameter given in (7) can be expressed as = (N − 1) , where
=(Cβ1− θ)T⊗ (Cβ1− θ)T
· vecCS∗−1CT−1/σ2.
Let ∂ /∂ vec(S∗)denote the p2-dimensional column vector whose ith component is the deriv-ative of with respect to the ith element of vec(S∗). It can be shown by applying the algebraic manipulation and matrix differentiation results that
∂ ∂vec(S∗)= S∗−1CTCS∗−1CT−1(Cβ1− θ) ⊗S∗−1CTCS∗−1CT−1(Cβ1− θ) /σ2.
For operational ease, the derivative is computed ignoring the symmetry of S∗. Then, it can be readily derived from the Cramer delta method that has the following large-sample distribution
∼·N (μ , ), (8) where μ = (Cβ1− θ)T C−1CT−1(Cβ1− θ)/σ2 and = (G ⊗ G)T(G⊗ G)/ (N− 1)σ4− μ2 /(N− 1)
with G= −1CT(C−1CT)−1(Cβ1− θ). For C = Ip and θ= 0p, it leads to the useful results that μ = βT1β1/σ2and = (β1⊗ β1)T(β1⊗ β1)/{(N − 1)σ4} − μ2 /(N− 1). Conse-quently, the conditional distribution of F∗| with = (N − 1) , and asymptotic distribution of described in the last two equations specify fully the proposed approximate distribution for test statistic F∗. It is clear under the null hypothesis that the distribution of F∗ remains as F (c, N− p − 1) under both fixed and random settings as in the special case of the multinormal distribution of Sampson (1974). Hence, the test is conducted by rejecting H0if F∗> Fc,N−p−1,α. However, the power function associated with the general linear hypothesis H0: Cβ1= θ versus H1: Cβ1= θ can be well approximated by synthesizing the results in (7) and (8) as
P{F∗> Fc,N−p−1,α}=. ∞ −∞P Fc, N− p − 1, (N − 1) > Fc,N−p−1,α · g( ) d , (9)
where g( ) is the normal pdf of defined in (8). The numerical computation of approximate power requires the evaluations of central and noncentral F cdfs and the one-dimensional inte-gration with respect to a normal pdf. Since all related functions are readily embedded in modern statistical packages such as the SAS system, no substantial computing efforts are required.
For the purpose of sample size determination, the approximate power function defined in (9) can be employed to calculate the sample size needed to test hypothesis H0: Cβ1= θ versus H1: Cβ1= θ in order to attain the specified power for the chosen significance level α, para-meter values β and σ2, and probability distribution f (X∗i). The necessary sample size can be found through a simple iterative search. To reduce the computational effort in the search process, the starting sample size can be selected from the following simplified examination. Consider the even stronger asymptotic results for S∗and SSE∗ that S∗converges in probability to and
SSE∗/{(N − p − 1)σ2} converges in probability to 1. It follows from the application of Slut-sky’s theorem that SSH∗/{(N − 1)σ2} converges in distribution to the chi-square distribution χ2(c, μ ), the noncentral chi-square distribution with c degrees of freedom and noncentrality parameter μ , where μ is defined in (8). More importantly, the distribution of the c· F∗ sta-tistic can be alternatively approximated by the distribution χ2(c, (N − 1)μ ). Therefore, the corresponding approximate power function is P{χ2(c, (N− 1)μ ) > χc,α2 }, where χc,α2 is the upper 100α percentage point of the central chi-square distribution χ2(c). Hence, the sample size, say NCS, required to achieve the specified power level is a one-time direct inversion of a non-central chi-square cdf. In general, the resulting sample size provides a close but smaller value than the desired outcome according to the proposed mixture of the noncentral F cdf in (9). Note that the probability P{F∗> Fc,N−p−1,α}, for fixed values of c, p, α and model parameters, is increasing in sample size N . Hence, by starting with sample size NCSfor N , it only requires a small number of incremental searches in order to find the minimum sample size that attains the nominal power.
It is noteworthy that the proposed approach avoids the need for a full specification of the joint distributional form of X∗i by only assuming the second- and fourth-order mixed central moments of the underlying distribution. However, the calculations are fairly straightforward for some well-known distributions and it may require more involved mathematical manipulations (integration or summation) for complex and nonstandard situations. On the contrary, the mean of X∗i is immaterial to the distribution of S∗and, more importantly, the suggested approximation. Additionally, it should be noted that the effect size δ given in (4) plays an important role in power and sample size determinations for fixed regression models. Owing to the proposed two-stage distribution approximation to the F∗statistic described in (7) and (8), there is no simple closed-form expression for the effect size in (9). However, it can be comprehended from the mean value μ of the random noncentrality parameter . Therefore, the computed value of μ is viewed as a pseudo effect size.
3. Numerical Examples
For illustrative purposes, we present in this section the power and sample size calculations for a random simple regression model and the moderated multiple regression or interaction re-gression model with two continuous predictor variables and their cross-product term.
First, the random simple regression (p= 1) of the form Y = β0+ Xβ1+ ε is investigated, where ε has a normal distribution N (0, σ2). Without loss of generality, both the intercept parame-ter β0and the variance σ2are taken to be 1. As suggested by a referee, we consider two classes of distributions for the explanatory variable X, namely standardized gamma and standardized Poisson distributions. Therefore, the mean and variance of X are identically μ= 0 and = 1, respectively. In order to investigate the finite-sample properties of the suggested procedure with respect to various shapes of distributions, the gamma distributions with shape parameter 9, 4 and 1 and scale parameter 1, denoted by gamma(9, 1), gamma(4, 1) and gamma(1, 1), and the Pois-son distributions with mean 9, 4 and 1, denoted by PoisPois-son(9), PoisPois-son(4) and PoisPois-son(1), are considered. Note that the skewness and kurtosis of the gamma(a, 1) distribution are 2/a1/2and 3+ 6/a, respectively. Hence, the actual values of skewness and kurtosis for the three prescribed gamma distributions are (0.67, 3.67), (1, 4.5) and (2, 9), respectively. In the case of Poisson distribution, the skewness and kurtosis of the Poisson(λ) distribution are 1/λ1/2and 3+ 1/λ, respectively. Thus, for the three Poisson distributions, the corresponding skewness and kurtosis are (0.33, 3.11), (0.5, 3.25) and (1, 4). For the test of slope coefficient (c= p = 1) H0: β1= 0 versus H1: β1= 0, it follows from (9) that the sample size needed to obtain the power 1− γ at the significance level α= 0.05 is the minimum number N such that the approximate power function P{F∗> F1,N−2,α}=. ∞ −∞P F1, N− 2, (N − 1) > F1,N−2,α · g( ) d ≥ 1 − γ, (10)
where ∼ N(μ. , ), μ = β12· /σ2and = β14· ( − 2)/[(N − 1)σ4]. For regression coefficient β1= 0.3, 0.4 and 0.5 and power level 1 − γ = 0.80, 0.90 and 0.95, the calculated sample sizes of the proposed method are presented in Tables1and2for the standardized gamma and standardized Poisson distributions of X, respectively. The results in the two tables reveal the general relation that sample sizes increase with increasing power and kurtosis, and decrease with increasing value of β1. To demonstrate the power computation, the precise achieved powers associated with the derived sample sizes are recalculated with the proposed approximation given in (10) for all cases. As expected, the resulting approximate powers are slightly larger than their corresponding nominal power levels. Specifically, the calculated sample sizes of the proposed method for gamma(9, 1) distribution with β1= 0.3 are 93, 124 and 152 for power 0.80, 0.90 and 0.95, respectively. The corresponding approximate powers are 0.8027, 0.9020 and 0.9500, and almost identical to 0.80, 0.90 and 0.95, respectively.
The second model under consideration is the simple interaction model: Y = β0+ Xβ1+ Zβ2+ XZβ3+ ε, where ε ∼ N(0, 1). The two predictors (X, Z) are jointly normally distributed with mean (0, 0), variance (1, 1) and correlation ρ. It is important to note that, although both X and Z are normally distributed, the interaction term XZ is obviously not a normal random variable. Therefore, the established methods for multinormal covariates are inappropriate for the power and sample size calculations of this interaction regression model. In this case, it can be shown that E[XZ] = ρ, E[X2Z] = E[XZ2] = 0 and V [XZ] = 1 + ρ2, see Aiken and West (1991, Appendix A). Therefore, we have
= E[H] = ⎡ ⎣ρ1 ρ1 00 0 0 1+ ρ2 ⎤ ⎦ , where H = ⎡ ⎣ X 2 XZ X(XZ− ρ) XZ Z2 Z(XZ− ρ) X(XZ− ρ) Z(XZ − ρ) (XZ − ρ)2 ⎤ ⎦ .
T ABLE 1. Calculated sample sizes, approximate po wers and simulated po wers of the p roposed method for random simple re gression models Y = β0 + Xβ 1 + ε with standardized g amma predictor (p = c = 1a n d α = 0 .05 ). Gamma (9 , 1 ) with kurtosis = 3 .67 Gamma (4 , 1 ) with kurtosis = 4 .5 G amma (1 , 1 ) with kurtosis = 9 N Approximate Simulated E rror N Approximate Simulated E rror N Approximate Simulated E rror po wer po w er po wer po w er po wer po w er (i) β1 = 0 .3 93 0.8027 0.8003 − 0 .0024 94 0.8039 0.8060 0 .0021 97 0.8001 0.7985 − 0 .0016 124 0.9020 0.9065 0 .0045 125 0.9017 0.9009 − 0 .0008 131 0.9012 0.9083 0 .0071 152 0.9500 0.9550 0 .0050 154 0.9506 0.9522 0 .0016 162 0.9500 0.9537 0 .0037 (ii) β1 = 0 .4 55 0.8058 0.8080 0 .0022 55 0.8006 0.7942 − 0 .0064 59 0.8012 0.8079 0 .0067 73 0.9036 0.8987 − 0 .0049 74 0.9031 0.9105 0 .0074 80 0.9013 0.9122 0 .0109 89 0.9503 0.9511 0 .0008 91 0.9512 0.9481 − 0 .0031 100 0.9510 0.9562 0 .0052 (iii) β1 = 0 .5 37 0.8050 0.8097 0 .0047 38 0.8079 0.8036 − 0 .0043 42 0.8083 0.8186 0 .0103 49 0.9032 0.8996 − 0 .0036 50 0.9022 0.9059 0 .0037 57 0.9037 0.9181 0 .0144 60 0.9510 0.9530 0 .0020 62 0.9521 0.9563 0 .0042 71 0.9505 0.9651 0 .0146 The v alues o f μ associated with the three distrib u tions are 0 .09, 0.16 and 0 .25 for the three coef ficient v alues β1 = 0 .3 , 0 .4 and 0.5, respecti v ely .
T ABLE 2. Calculated sample sizes, approximate po wers and simulated po wers of the p roposed method for random simple re gression models Y = β0 + Xβ 1 + ε with standardized Poisson p redictor (p = c = 1a n d α = 0 .05 ). Poisson(9) with kurtosis = 3 .11 Poisson(4) with kurtosis = 3 .25 Poisson(1) with kurtosis = 4 N Approximate Simulated E rror N Approximate Simulated E rror N Approximate Simulated E rror po wer po w er po wer po w er po wer po w er (i) β1 = 0 .3 92 0.8004 0.7986 − 0 .0018 93 0.8042 0.7978 − 0 .0064 93 0.8015 0.7952 − 0 .0063 123 0.9014 0.9051 0 .0037 123 0.9010 0.8987 − 0 .0023 124 0.9010 0.9036 0 .0026 151 0.9500 0.9481 − 0 .0019 152 0.9509 0.9494 − 0 .0015 153 0.9505 0.9500 − 0 .0005 (ii) β1 = 0 .4 54 0.8019 0.8001 − 0 .0018 54 0.8010 0.8021 0 .0011 55 0.8037 0.8032 − 0 .0005 72 0.9027 0.9038 0 .0011 72 0.9019 0.9050 0 .0031 73 0.9019 0.9000 − 0 .0019 88 0.9503 0.9506 0 .0003 89 0.9519 0.9526 0 .0007 90 0.9511 0.9526 0 .0015 (iii) β1 = 0 .5 37 0.8103 0.8051 − 0 .0052 37 0.8090 0.8098 0 .0008 37 0.8019 0.8040 0 .0021 48 0.9018 0.9033 0 .0015 48 0.9006 0.9012 0 .0006 49 0.9005 0.9080 0 .0075 59 0.9511 0.9510 − 0 .0001 59 0.9503 0.9495 − 0 .0008 61 0.9520 0.9547 0 .0027 The v alues o f μ associated with the three distrib u tions are 0 .09, 0.16 and 0 .25 for the three coef ficient v alues β1 = 0 .3 , 0.4 and 0.5, respecti v ely .
Also, can be expressed as = ⎡ ⎣12 24 35 3 5 6 ⎤ ⎦ = E[H ⊗ H] = E ⎡ ⎣ X 2H XZH X(XZ− ρ)H XZH Z2H Z(XZ− ρ)H X(XZ− ρ)H Z(XZ − ρ)H (XZ − ρ)2H ⎤ ⎦ , where 1= ⎡ ⎣3ρ3 1+ 2ρ3ρ 2 00 0 0 3+ 7ρ2 ⎤ ⎦ , 2= ⎡ ⎣ 3ρ 1+ 2ρ 2 0 1+ 2ρ2 3ρ 0 0 0 7ρ+ 3ρ3 ⎤ ⎦ , 3= ⎡ ⎣ 0 0 3+ 7ρ 2 0 0 7ρ+ 3ρ3 3+ 7ρ2 7ρ+ 3ρ3 0 ⎤ ⎦ , 4= ⎡ ⎣1+ 2ρ 2 3ρ 0 3ρ 3 0 0 0 3+ 7ρ2 ⎤ ⎦ , 5= ⎡ ⎣ 0 0 7ρ+ 3ρ 3 0 0 3+ 7ρ2 7ρ+ 3ρ3 3+ 7ρ2 0 ⎤ ⎦ and 6= ⎡ ⎣ 3+ 7ρ 2 7ρ+ 3ρ3 0 7ρ+ 3ρ3 3+ 7ρ2 0 0 0 9+ 42ρ2+ 9ρ4 ⎤ ⎦ .
The evaluations of for this interaction model are more involved than those in the second model. The conditional distribution properties of Z given X and high-order moments of a standard nor-mal distribution (E[X6] = 15 and E[X8] = 105) are required to carry out the calculations. For illustration, the coefficient parameters are set as (β0, β1, β2, β3)= (1, 0.1, 0.3, 0.25). Two hy-pothesis tests are investigated in this numerical demonstration. They are the tests of overall ef-fects (c= 3) and interaction effect (c = 1) with the null hypotheses H0: β1= β2= β3= 0 and H0: β3= 0, respectively. In a similar fashion, the proposed approach is employed to perform the sample size and the corresponding approximate power calculations for testing the specified hy-pothesis with significance level α= 0.05 and nominal power (0.80, 0.90, 0.95). These numerical results are presented in Table3for three different values of ρ= 0.3, 0.5 and 0.7.
It is important to note that the major analytical justification considered here applies large-sample approximation to the distribution of the F∗statistic. In order to assess the finite-sample accuracy of the proposed approach, simulation studies are conducted next. With given sample size and model configuration, an estimate of the true power or simulated power is then com-puted through simulation of 10,000 replicate data sets. For each replicate, N sets of explanatory variables are generated from the selected distribution. These values in turn determine the mean responses for generating N normal outcomes with the underlying regression model. Then the test statistic is computed and the simulated power is the proportion of the 10,000 replicates whose F∗test statistic values exceed the critical value Fc,N−p−1,α. The adequacy of the proposed sam-ple size formula is determined by the difference (simulated power–approximate power) between the simulated power and approximate power specified above. All calculations are performed us-ing programs written with SAS/IML (SAS Institute,2003). Detailed numerical results of the simulation studies are reported in Tables1–2 and3 for the two models, respectively. For the simple regression models, under the standardized gamma(9, 1) predictor variable situation with β1= 0.3, the simulated powers are 0.8003, 0.9065 and 0.9550 for the three different power levels= 0.80, 0.90 and 0.95, respectively. Thus, the differences or errors between simulated
TABLE3.
Calculated sample sizes, approximate powers and simulated powers of the proposed method for random multiple regres-sion model Y= β0+ Xβ1+ Zβ2+ XZβ3+ ε (p = 3, α = 0.05) with standard normal error.
Test of overall effects (c= 3) Test of interaction effect (c= 1)
N Approximate Simulated Error N Approximate Simulated Error
power power power power
(i) ρ= 0.3 70 0.8029 0.8067 0.0038 127 0.8013 0.7970 −0.0043 91 0.9002 0.9032 0.0030 171 0.9010 0.8985 −0.0025 111 0.9502 0.9558 0.0056 212 0.9503 0.9494 −0.0009 (ii) ρ= 0.5 65 0.8049 0.8015 −0.0034 114 0.8012 0.7930 −0.0082 85 0.9017 0.9103 0.0086 154 0.9007 0.8933 −0.0074 104 0.9508 0.9557 0.0049 192 0.9505 0.9515 0.0010 (iii) ρ= 0.7 60 0.8067 0.8148 0.0081 99 0.8010 0.7982 −0.0028 79 0.9028 0.9144 0.0116 135 0.9015 0.9010 −0.0005 97 0.9511 0.9583 0.0072 169 0.9510 0.9535 0.0025
The values of μ associated with the two tests are (0.1861, 0.0681), (0.2081, 0.0781) and (0.2351, 0.0931) for the three correlation values ρ= 0.3, 0.5, and 0.7, respectively.
powers and approximate powers are 0.8003− 0.8027 = −0.0024, 0.9065 − 0.9020 = 0.0045 and 0.9550− 0.9500 = 0.0050, respectively. Similarly, all other results in Tables1–3 are ob-tained.
Examination of Table 1 shows that the absolute errors associated with the standardized gamma(9, 1) and gamma(4, 1) predictor distributions do not exceed 0.01 for all combinations of β1and power levels. However, the results for the standardized gamma(1, 1) distribution vary with the value of β1and power levels. Specifically, the cases associated with β1= 0.3 are less sensitive to the influence of the outsized skewness 2 and kurtosis 9 of the standardized gamma(1, 1) distri-bution than those for β1= 0.4 and β1= 0.5. Obviously, the errors 0.0103, 0.0144 and 0.0146 for the three power levels of β1= 0.5 are greater than 0.01. Hence, the accuracy of the large-sample power approximation is not as satisfactory as other circumstances for strongly skewed gamma distributions, especially for small samples. Nonetheless, the relative performance of the proposed method for Poisson distributions in Table2is excellent even for the most skewed Poisson(1) dis-tribution with comparatively small sample sizes.
With respect to the second model, the results in Table3suggest that there is a close agree-ment between the simulated power and the approximate power because the absolute errors are less than 0.01. The only exception is 0.0116 which is associated with the test of overall effects for correlation ρ= 0.7. Overall, the accuracy of the proposed approach increases slightly with the sample size, and varies marginally with the model configurations. According to these find-ings, the performance of the proposed method appears to be excellent for the range of random regression specifications considered here. It is important to note that, in the context of moderated multiple regression, the test for the existence of moderator effect is examined by the significance of regression coefficient β3of the cross-product term in the simple interaction model. As noted in the review by Aguinis et al. (2005), however, it has been widely recognized that moderated multiple regression analyses have suffered low statistical power in detecting moderator effects. Therefore, the proposed power formula can be employed to determine the minimum sample size required for testing the hypothesis H0: β3= 0 with specified model configurations, significance level, and nominal power.
T ABLE 4. Calculated sample sizes, approximate po wers and simulated po wers of the p roposed method for random multiple re gression model Y = β0 + Xβ 1 + Zβ 2 + XZ β3 + ε (p = 3, α = 0 .05 ) with standardized uniform error . T est of o v erall ef fects (c = 3) T est of interaction ef fect (c = 1) N Simulated E rror A pproximate Simulated E rror N Simulated E rror A pproximate Simulated E rror α po wer po w er α po wer po w er (i) ρ = 0 .3 70 0.0472 − 0 .0028 0.8029 0.7999 − 0 .0030 127 0.0490 − 0 .0010 0.8013 0.7833 − 0 .0180 91 0.0500 0 .0000 0.9002 0.9047 0 .0045 171 0.0484 − 0 .0016 0.9010 0.8903 − 0 .0107 111 0.0468 − 0 .0032 0.9502 0.9554 0 .0052 212 0.0472 − 0 .0028 0.9503 0.9523 0 .0020 (ii) ρ = 0 .5 65 0.0457 − 0 .0043 0.8049 0.8139 0 .0090 114 0.0471 − 0 .0029 0.8012 0.7887 − 0 .0125 85 0.0508 0 .0008 0.9017 0.9110 0 .0093 154 0.0426 − 0 .0064 0.9007 0.9044 0 .0037 104 0.0459 − 0 .0041 0.9508 0.9612 0 .0104 192 0.0484 − 0 .0016 0.9505 0.9469 − 0 .0036 (iii) ρ = 0 .7 60 0.0499 − 0 .0001 0.8067 0.8123 0 .0056 99 0.0504 0 .0004 0.8010 0.7912 − 0 .0098 79 0.0484 − 0 .0016 0.9028 0.9155 0 .0127 135 0.0501 0 .0001 0.9015 0.8994 − 0 .0021 97 0.0476 − 0 .0024 0.9511 0.9638 0 .0127 169 0.0510 0 .0010 0.9510 0.9516 0 .0006
T ABLE 5. Calculated sample sizes, approximate po wers and simulated po wers of the p roposed method for random multiple re gression model Y = β0 + Xβ 1 + Zβ 2 + XZ β3 + ε (p = 3, α = 0 .05 ) with standardized g amma error . T est of o v erall ef fects (c = 3) T est of interaction ef fect (c = 1) N Simulated E rror A pproximate Simulated E rror N Simulated E rror A pproximate Simulated E rror α po wer po w er α po wer po w er (i) ρ = 0 .3 70 0.0507 0 .0007 0.8029 0.8115 0.0086 127 0.0509 0 .0009 0.8013 0.7936 − 0 .0077 91 0.0510 0 .0010 0.9002 0.9062 0.0060 171 0.0480 − 0 .0020 0.9010 0.8958 − 0 .0052 111 0.0501 0 .0001 0.9502 0.9551 0.0049 212 0.0497 − 0 .0003 0.9503 0.9512 0 .0009 (ii) ρ = 0 .5 65 0.0514 0 .0014 0.8049 0.8051 0.0002 114 0.0511 0 .0011 0.8012 0.7862 − 0 .0150 85 0.0490 − 0 .0010 0.9017 0.9051 0.0034 154 0.0478 − 0 .0022 0.9007 0.9013 0 .0006 104 0.0522 0 .0022 0.9508 0.9573 0.0065 192 0.0489 − 0 .0011 0.9505 0.9556 0 .0051 (iii) ρ = 0 .7 60 0.0503 0 .0003 0.8067 0.8122 0.0055 99 0.0495 − 0 .0005 0.8010 0.7900 − 0 .0110 79 0.0511 0 .0011 0.9028 0.9140 0.0112 135 0.0485 − 0 .0015 0.9015 0.9048 0 .0033 97 0.0502 0 .0002 0.9511 0.9625 0.0114 169 0.0501 0 .0001 0.9510 0.9559 0 .0049
As pointed out by a referee, Anderson (1999) showed that the coefficients estimator within the multivariate multiple regression framework has an asymptotic multinormal distribution when the errors and predictors are mutually independently distributed, irrespective of whether they are normal. Since the multiple regression model considered here is a special case of the multi-variate multiple regression, the asymptotic normality property of the multiple regression coeffi-cient estimator can readily be established. However, Anderson (1999) did not explicitly discuss the distribution of the associated F statistic and power calculation. It is important to note that the conditional distribution of the F∗ statistic given in (6) is no longer necessarily an exact F distribution. Nonetheless, the proposed large-sample approximation for the conditional normal regression model given in (7–9) can be applied to the situation of nonnormal errors and pre-dictors. To examine the robust issues of the proposed method against the extra complication of nonnormal errors, we have conducted a numerical evaluation for the simple interaction model: Y = β0+ Xβ1+ Zβ2+ XZβ3+ ε, where ε has a standardized uniform(0, 1) or standardized gamma(5, 1) distribution. For ease of exposition, the parameter settings are the same as those in Table3, and the results are presented in Tables4and5. The simulated Type I error rate and power are compared with the nominal α= 0.05 and power level, respectively. Notably, all the absolute errors between simulated α and nominal value 0.05 are less than 0.01. In addition, the discrepancies between simulated powers and approximate powers calculated with the proposed power function (9) are slightly larger than those in Table 3with standard normal errors. How-ever, the performance seems completely acceptable, given the many unknowns in study planning. Therefore, the suggested procedures for conditional normal regression models with arbitrary dis-tribution configurations for explanatory variables are not seriously affected by mild departures from the normality assumption of errors.
4. Conclusions
Procedures for power and sample size determinations in fixed regression models have been developed for years but none seems to have provided a comprehensive treatment or guideline for the calculations of power and sample sizes in the framework of random regression mod-els. Within the context of random regression models, the current results are mainly under the situation of the normality assumption for explanatory variables. A natural generalization to in-corporate other distributions of explanatory variables is essential to researchers for performing power and sample size calculations in practice. This paper discusses a feasible solution to this issue by providing both theoretical justification and numerical examination for the proposed uni-fied approach. With this direct extension, one can perform power and sample size calculations in multiple regression models with any discrete and/or continuous distributions of explanatory variables. The remarkable performance for power and sample size calculations reveals that the proposed method may find useful applications in subsequent random regression analysis. No-tably, the suggested methodology is applicable for the prominent moderated multiple regression models containing a continuous predictor and moderator.
References
Aguinis, H., Beaty, J.C., Boik, R.J., & Pierce, C.A. (2005). Effect size and power in assessing moderating effects of categorical variables using multiple regression: A 30-year review. Journal of Applied Psychology, 90, 94–107. Aiken, L.S., & West, S.G. (1991). Multiple regression: Testing and interpreting interactions. Thousand Oaks: Sage. Algina, J., & Olejnik, S. (2000). Determining sample size for accurate estimation of the squared multiple correlation
coefficient. Multivariate Behavioral Research, 35, 119–137.
Anderson, T.W. (1999). Asymptotic distribution of the reduced rank regression estimator under general conditions. An-nals of Statistics, 27, 1141–1154.
Gatsonis, C., & Sampson, A.R. (1989). Multiple correlation: Exact power and sample size calculations. Psychological Bulletin, 106, 516–524.
Kelley, K., & Maxwell, S.E. (2003). Sample size for multiple regression: Obtaining regression coefficients that are accurate, not simply significant. Psychological Methods, 8, 305–321.
Mendoza, J.L., & Stafford, K.L. (2001). Confidence interval, power calculation, and sample size estimation for the squared multiple correlation coefficient under the fixed and random regression models: A computer program and useful standard tables. Educational and Psychological Measurement, 61, 650–667.
Muirhead, R.J. (1982). Aspects of multivariate statistical theory. New York: Wiley.
Raudenbush, S.W., & Liu, X. (2000). Statistical power and optimal design for multisite randomized trials. Psychological Methods, 5, 199–213.
Raudenbush, S.W., & Liu, X. (2001). Effects of study duration, frequency of observation, and sample size on power in studies of group differences in polynomial change. Psychological Methods, 6, 387–401.
Rencher, A.C. (2000). Linear models in statistics. New York: Wiley.
Sampson, A.R. (1974). A tale of two regressions. Journal of the American Statistical Association, 69, 682–689. SAS Institute (2003). SAS/IML user’s guide, Version 8. Cary, NC: SAS Institute Inc.
Shieh, G. (2003). A comparative study of power and sample size calculations for multivariate general linear models. Multivariate Behavioral Research, 38, 285–307.
Shieh, G. (2005). Power and sample size calculations for multivariate linear models with random explanatory variables. Psychometrika, 70, 347–358.
Shieh, G. (2006). Exact interval estimation, power calculation and sample size determination in normal correlation analy-sis. Psychometrika, 71, 529–540.
Steiger, J.H., & Fouladi, R.T. (1992). R2: A computer program for interval estimation, power calculations, sample size es-timation, and hypothesis testing in multiple regression. Behavioral Research Methods, Instruments, and Computers, 24, 581–582.
Manuscript received 1 SEP 2005 Final version received 18 DEC 2006 Published Online Date: 7 JUL 2007
