Introduction

Interactive analysis has received fairly extensive attention in biological and medical sciences. The occurrence of interaction implies that the influence of a predictor variable on the response variable varies as a function of the level of a treatment variable. In general, the existence and magnitude of interaction effects can be measured and appraised by the statistical methods for analyzing the product terms of treatment and predictor variables. Alternatively, the procedure for evaluating the interaction effects between treatment and predictor variables is methodologically identical to that for testing the equality of regression slope coefficients in two or more regression lines. The common testing procedures of multiple linear regression described by Fleiss [1], Kleinbaum et al. [2], and Kutner et al. [3] can be readily applied for the detection of different forms of interaction effects. Among various pedagogical development and exposition of interaction studies, particular emphasis has been paid to the fundamental attributes of two treatment groups, such as Hewett and Lababidi [4], Hill and Padmanabhan [5], Spurrier, Hewett, and Lababidi [6], and Tsutakawa and Hewett [7], among others.

The traditional multiple linear regression adopts the fundamental assumptions of independence, normality, and constant variance. Accordingly, homogeneity of error variance is one of the essential assumptions for the appropriate use of common test procedures. The detrimental effects of heterogeneous error variance on the regular test procedures have been investigated in Alexander and DeShon [8], DeShon and Alexander [9], Dretzke, Levin, and Serlin [10], Overton [11], and Shieh [12]. The numerical findings indicated that the resulting Type I error rate and power levels may be considerably affected when group sample sizes are equal, and severely distorted when group sample sizes are unequal. Therefore, the traditional tests of interaction effects are not robust to the heterogeneity of error variance for evaluating regression slope differences across groups. It is emphasized in Dretzke et al. [10] and Shieh [12] that the importance of choosing a proper test procedure for two regression slope differences when group error variances are unequal. Moreover, under the prevalent scenario of variance heterogeneity, a theory-based approach to power calculation for assessing the existence of treatment by predictor interactions remains unavailable.

To enhance the adoption of appropriate techniques for statistical inference and research design, this article aims to present power and sample size procedures for the well-recognized Welch [13] test of parallelism of two regression lines with heterogeneous variances. Accordingly, the suggested techniques are direct extensions of those considered in Shieh [14] under a homogeneous variance setting, just as the Welch [13] procedure to the traditional two-sample t test for mean comparisons. Four important aspects of this study should be pointed out. First, the primitive Welch solution to the Behrens-Fisher problem for comparing the means of two groups under heterogeneity of variance has been widely discussed in the literature (Kim and Cohen [15]). Although there exist several viable alternatives to mitigate the impact of violating the assumption of homogeneous error variance, the accurate performance and computational ease demonstrated in Dretzke et al. [10] and Shieh [12] ensure that the extended Welch procedure for regression slope equality can be recommended. Second, the actual values of predictors cannot be known in advance before conducting a research study just as the primary responses. In order to account for the stochastic nature of predictors, two different general or unconditional distributions of the extended Welch statistic are described to obtain feasible power functions and sample size formulas. Third, comprehensive numerical assessments were conducted to illustrate the discrepancy between two potential approaches under a wide range of model configurations. The empirical findings revealed that the more elaborate procedure generally outperforms the alternative simplified method for correctly accommodating the predictor features. Fourth, to ease the computational demands of the proposed power and sample size calculations, both SAS and R computer algorithms are presented. The program codes are computationally efficient and will give accurate outputs when all the necessary information is properly specified. With the most plausible model configurations, the recommended techniques are practically useful for power analysis and sample size determination in planning interaction studies.

Methods

Consider the two simple linear regression models with unknown and possibly unequal coefficient parameters and error variances: (1) where ε_1j and ε_1k are iid and ) random variables, respectively, j = 1, …, N₁, and k = 1, …, N₂. Note that Shieh [14] focused on the scenario with homogeneous error variance . Moreover, the regression model with unequal slopes in Eq 1 can be formulated as an interactive model using a dichotomous treatment variable M: (2) Evidently, the slope difference between the two regression lines β_1D = β₁₁−β₁₂ is also the coefficient associated with the cross-product term between the predictor and treatment variables. The existence and magnitude of β_1D indicates the influences of a dichotomous treatment indicator on the direction and/or the strength of the relationship between the response and predictor variables.

For the purpose of detecting the interaction effect, the problem is naturally concerned with the hypothesis (3) and the distributional properties of the corresponding least squares estimators of slope coefficients and . Under the model formulation in Eq 1, the standard results show that and where , and and are the respective sample means of the X_1j and X_2k observations. These results lead to where . Also, and are the usual unbiased estimators of and where SSE₁ and SSE₂ are the separate group error sum of squares, c₁ = N₁−2, and c₂ = N₂−2. Moreover, and , where χ²(c₁) and χ²(c₂) are chi-square distribution with c₁ and c₂ degrees of freedom, respectively.

To detect the non-parallelism of two regression lines or the interactive effects of a dichotomous treatment variable between the predictor and response variables in terms of H₀: β_1D = 0 versus H₁: β_1D ≠ 0, Welch [13] proposed an extension of the well-known approximate degrees of freedom solution to the Behren-Fisher problem of comparing the difference between two means. Accordingly, the extended Welch statistic has the form: (4) Under the null hypothesis H₀: β_1D = 0, Welch [13] considered the approximate distribution for T*: (5) where is a t distribution with degrees of freedom and The null hypothesis is rejected at the significance level α if (6) where is the 100(1−α/2) percentile of the t distribution with degrees of freedom .

The statistical inferences about the interactive effect are based on the conditional distribution of the continuous outcomes {X_1j, j = 1, …, N₁} and {X_2k, k = 1, …, N₂}. Consequently, the particular results would depend on the observed values of the predictor variables. At the design stage, however, the actual values of both the response and predictor variables cannot be determined. Hence, it is more appropriate to consider the random or unconditional framework as demonstrated in Binkley and Abbot [16], Cramer and Appelbaum [17], Gatsonis and Sampson [18], and Sampson [19]. Although the unconditional properties of the test procedure are relatively more involved, the associated hypothesis testing and parameter estimation remain the same under both conditional and unconditional modeling settings. The fundamental distinction between the two modeling scenarios is essential when performing power and sample size calculations. It is of theoretical importance to recognize the random characteristic of the predictor variables and to appraise the unconditional property of the test statistic over the possible values of the predictors.

In order to explicate the critical notion of accommodating the distributional properties of the predictor variables, two different power functions are presented in Appendix. The simplified t method only utilizes the partial information of predictor variances. The associated power function Ψ_ST defined in Eq A3 has the advantages of statistical and computational simplifications. On the other hand, the mixed t approach considers the improved power function Ψ_MT given in Eq A6 and accommodates the full distributional features of the predictors. Despite the close resemblance between the two power formulas, the corresponding behavior may differ substantially for finite samples. In practice, a research study requires adequate statistical power and sufficient sample size to detect scientifically credible effects. For the purpose of planning interaction design, the prescribed two power formulas can be employed to determine the sample sizes N₁ and N₂ needed to attain the specified power through a simple iterative search for the chosen significance level and model configurations. Accordingly, it is of practical interest and theoretical importance to examine the underlying characteristics of the two approaches in power and sample size calculations.

Simulation study

Due to the complexity of the power and sample size problem, a complete theoretical evaluation of the associated issues is not feasible. To clarify the similarities and differences of the prescribed two different distributions for the extended Welch statistic, extensive numerical appraisals were conducted for their performance in power and sample size calculations. The investigation was carried out in two steps. The first step involved sample size determinations across a wide range of manipulated model configurations. In the second step, a Monte Carlo simulation study was performed to compare the accuracy of the alternative sample size formulas through power assessments under the design characteristics specified in the first step.

Sample size determinations

The determination of necessary sample size for testing regression slope differences requires detailed specifications of the magnitudes of slope coefficients, error variances, predictor variances, nominal power, and significance level. Notably, the impact of each of these critical factors on sample size calculations not only differs but also varies with the coexisting magnitudes of other components. To provide a systematic evaluation, the empirical examinations are administered by fixing all but one of the following critical attributes and by varying a single component in the appraisal: (1) slope differences, (2) sample size ratios, (3) predictor variances, and (4) error variances. Specifically, the magnitudes of slope coefficients are chosen as β₁₁ = {0.50, 0.75, 1.00} and β₁₂ = 0 to represent small, medium, and large slope differences with β_1D = {0.50, 0.75, 1.00}. Both balanced and unbalanced designs are considered with sample size ratio r = N₂/N₁ = 1 and 3. The predictor variances are assumed to have three different structures: = {1, 1}, {1, 4}, and {4, 1}. Moreover, a total of five patterns of error variances are evaluated: = {1, 1}, {1, 2}, {1, 4}, {2, 1}, and {4, 1}. The diverse settings not only include both balanced and unbalanced schemes, but also generate direct- and inverse-pairing with variance component patterns. Note that the empirical investigation involves a total of 90 combined conditions. It should be emphasized that these model configurations are designated to represent the possible extent of characteristics in practical study without unrealistic or excessive settings. In fact, similar setups were considered in Overton [11] and Shieh [12]. Throughout these numerical comparisons, the significance level and nominal power are set as α = 0.05 and 1−β = 0.80, respectively.

With the aforementioned specifications, the necessary sample sizes can be computed by the supplemental SAS/IML (SAS Institute [20]) or R (R Development Core Team [21]) programs for the simplified t and mixed t methods with the power functions defined in Eqs A3 and A6, respectively. The resulting sample sizes {N₁, N₂} associated with the slope difference β_1D = {0.05, 0.75, 1.00} in combination with sample size ratio r = {1, 3} are summarized in Tables 1–6, respectively. Accordingly, there are 15 combined cases of and in each table. For ease of illustration, the total sample size N_T and the working effect size δ* given in Eq A4 are also presented. As expected, the computed sample sizes increase with smaller magnitude of slope difference β_1D or effect size δ*. Also, the necessary sample sizes associated with larger error variances are consistently larger than those of the smaller error variances when all other factors are held constant. However, the situation is reversed in terms of the predictor variances. In other words, the scenarios with larger incurred smaller sample sizes than those with smaller predictor variances when all other characteristics remain the same.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Computed sample size, estimated power, and simulated power when the slope difference β_1D = 0.50, sample size ratio r = 1, Type I error α = 0.05, and nominal power 1−β = 0.80.

https://doi.org/10.1371/journal.pone.0207745.t001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Computed sample size, estimated power, and simulated power when the slope difference β_1D = 0.50, sample size ratio r = 3, Type I error α = 0.05, and nominal power 1−β = 0.80.

https://doi.org/10.1371/journal.pone.0207745.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Computed sample size, estimated power, and simulated power when the slope difference β_1D = 0.75, sample size ratio r = 1, Type I error α = 0.05, and nominal power 1−β = 0.80.

https://doi.org/10.1371/journal.pone.0207745.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Computed sample size, estimated power, and simulated power when the slope difference β_1D = 0.75, sample size ratio r = 3, Type I error α = 0.05, and nominal power 1−β = 0.80.

https://doi.org/10.1371/journal.pone.0207745.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Computed sample size, estimated power, and simulated power when the slope difference β_1D = 1.00, sample size ratio r = 1, Type I error α = 0.05, and nominal power 1−β = 0.80.

https://doi.org/10.1371/journal.pone.0207745.t005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Computed sample size, estimated power, and simulated power when the slope difference β_1D = 1.00, sample size ratio r = 3, Type I error α = 0.05, and nominal power 1−β = 0.80.

https://doi.org/10.1371/journal.pone.0207745.t006

Also, the unbalanced designs often demand larger sample size to attain the same nominal power than the balanced designs. But there are some exceptions, for example, Cases 3, 11, 12, and 13 in Tables 1 and 2 signify the notable situations. Specifically, the total sample size N_T = 320 and 324 are slightly larger than those of N_T = 300 and 304 for the two methods associated with Case 3 in Tables 1 and 2, respectively. More importantly, the sample sizes of the simplified t method are relatively smaller than those of the mixed t approach. The only one situation with the same sample sizes is associated with {N₁, N₂} = {24, 72} of Case 13 in Table 4 when β_1D = 0.75, = {4, 1}, and = {1, 4}. However, with the two different power functions Ψ_ST and Ψ_MT, the corresponding estimated powers are 0.8164 and 0.8028, respectively. Note that if problems or deficiencies were to be seen with the alternative power and sample size procedures, they would be most apparent with small sample sizes. Thus, some of these sample sizes are designated to be smaller than those in related studies. In order to evaluate the accuracy of the two power functions, the estimated power or computed power are also listed. Because of the underlying metric of integer sample sizes, the attained values are slightly larger than the nominal level for both procedures.

Results

Conclusions and discussion

Assessing interaction effects of categorical treatment variables or testing the equality of regression slopes is a fundamental procedure in biology, medicine, and other fields of science. The present study focuses on the two-group case because it represents the great majority of interactive studies. In view of the frequent violation of the homogeneity of error variance assumption in applications, the extended Welch approach has been proposed as a vital alternative to the ordinary least squares method for the detection of interaction effect between a dichotomous treatment variable and a continuous predictor. However, the associated power computation and sample size determination for interaction analysis remain unavailable. To alleviate the difficulty in detecting interaction due to insufficient statistical power, this article seeks to contribute to the literature by presenting practical solutions to power and sample size calculations for the extended Welch tests of regression slope difference under variance heterogeneity. Accordingly, two distinct techniques are described to emphasize and account for the stochastic nature of the predictor variables for approximating the general distribution and power behavior of the extended Welch statistic. In addition to the detailed theoretical explications for the inherent relationship and distinct property of the two approaches, numerical investigations are conducted to reveal the relative performance of the contending procedures as a reliable technique to accommodate the predictor features that make interaction research particularly distinctive. The suggested approach has the important advantages over the simplified method in theoretical justification, great accuracy, and wide applicability. Due to the limitation of existing software packages and the practical usefulness of the proposed methodology, computer algorithms are presented to implement the recommended power calculations and sample size determinations.

Within the context of statistical point estimation, unbiasedness is an essential principle in the selection of desirable point estimators. However, the unbiased criterion does not guarantee the chance of obtaining a correct estimate. For example, it is well known that the sample mean is unbiased of the population mean. Although the probability of making a correct estimate is zero for the case of continuous variables, the sample mean was not deprived of its role as a reliable and established estimator. For the problem of sample size determination, some researchers may suggest that it makes no difference to apply a simple formula given the many unknowns in the planning stage of research design. Nevertheless, the underlying principle of research design planning is to conduct advance analysis based on the best knowledge of the model configurations. It is prudent to adopt a feasible procedure that incorporates into power and sample size calculations as much information about the model configurations as possible. Note that hierarchical linear models (Raudenbush & Bryk [26]) have the distinct structure for accommodating the notion of random coefficients in the traditional linear models. Here, the stochastic nature of predictor variables is taken into account for the described power and sample size methodology for interaction research. Essentially, the developed algorithms will yield accurate power appraisal and sample size determination provided that all the required features are suitably specified.

In addition to the relaxed scenario of variance heterogeneity, it should be stressed that the recommended power and sample size approaches also require the independence and normality assumptions. The robust feature clearly depends on the extent of how poorly the underlying response and predictor distributions diverge from the normality setting. However, Poncet et al. [27] and Schmider et al. [28] reinvestigated the robustness of two-sample t test and ANOVA F test against violations of the normal distribution assumption by Monte Carlo simulations. Their results gave strong support for the robustness of the t and F tests under application of non-normally distributed data. Although they did not evaluate the Welch method, the updated findings in some sense reinforce the documented robustness for the Welch-type procedures under mild departures from normality assumption. More importantly, Lix and Keselman [29] noted that the Welch-type tests preserve reasonably good Type I error control and compare favorably with the traditional method in terms of statistical power for many types and degrees of non-normality likely to be encountered in practical research. Moreover, the established class of generalized linear models offers an excellent and distinct method for analyzing non-normal data. Related details and follow-up procedures can be found in McCullagh and Nelder [30] and McCulloch, Searle, and Neuhaus [31].

Supporting information