Logistic regression

A logistic regression model describes the relationship between a binary response variable Y and one or more independent variables (predictors) X_i. The variables X_i can be discrete or continuous. In a simple logistic regression with one predictor X the assumption is that the event rate P = prob(Y = 1) depends on X in the following way:

P = e^{β0 +β1 X}/(1 + e^{β0 +β1 X})

For β₁ ≠ 0 and a continuous variable X this formula describes a smooth S-shaped transition of the probability for Y = 1 from 0 to 1 (β₁ > 0) or from 1 to 0 (β₁ < 0) with increasing X. This transition gets steeper with increasing β₁. Rearranging the formula leads to: log(P/(1 − P)) = β₀ + β₁ X. This shows that the logarithm of the odds P/(1 − P), also called a logit, on the left side of the equation is linear in X. Here, β₁ is the slope of this linear relationship.

The interesting question is whether predictor X_i is related to Y or not. Thus, in a simple logistic regression model , the null and alternative hypothesis for a two-sided test are:

H₀ : β₁ = 0,
H₁ : β₁ ≠ 0.

The procedures implemented in G*Power for this case estimates the power of the Wald test. The standard normally distributed test statistic of the Wald test is:
where βˆ₁ is the maximum likelihood estimator for parameter β₁ and var(βˆ₁) the variance of this estimate.

In a multiple logistic regression model with log(P/(1 − P)) = β₀ + β₁ X₁ + · · · + β_p X_p we test the effect of a specific predictor in the presence of other predictors. In this case the null hypothesis is

H₀ : [β₁, β₂, . . . , β_p] = [0, β₂ , . . . , β_p]

and the alternative hypothesis is

H₁ : [β₁, β₂, . . . , β_p] = [βˆ, β2 , . . . , βp ], where βˆ ≠ 0.

Effect size index

In the simple logistic model the effect of X on Y is given by the size of the parameter β₁. Let p₁ denote the probability of an event under H₀, that is exp(β₀) = p₁/(1 − p₁)- Let p₂ denote the probability of an event under H₁ at X = 1, that is exp(β₀ + β₁) = p₂/(1 − p₂).

Then exp(β₀ + β₁)/exp(β₀) = exp(β₁) = [p₂/(1 − p₂)]/[p₁/(1 − p₁)] := odds ratio OR, which implies β1 = log[OR].

Given the probability p₁ the effect size is specified by the odds ratio (OR) or directly by p₂. Click Options in the main window, and then select Input effect size as... either Odds ratio or Two probabilities.

Selecting Odds ratio will reveal the Odds ratio input field in the main window. Selecting Two probabilities will reveal the Pr(Y=1|X=1) H1 input field in the main window.

The probability p₁ is specified in the Pr(Y=1|X=1) H0 input field in the main window. OR is specified in the Odds ratio input field and, alternatively, p₂ is specified in the Pr(Y=1|X=1) H1 input field.

Setting p₂ = p₁ or equivalently OR = 1 implies β₁ = 0 and thus an effect size of zero. An effect size of zero must not be used in a priori analyses. Apart from these values, the following additional inputs are needed.

R2 other X.

In models with more than one covariate, the influence of the other covariates X₂, . . . , X_p on the power of the test can be taken into account by using a correction factor. This factor depends on the proportion R² = ρ²_1·23...p of the variance of X₁ explained by the regression relationship with X₂, . . . , X_p. If N is the sample size considering X₁ alone, then the sample size in a setting with additional covariates is: N' = N/(1 − R²). This correction for the influence of other covariates has been proposed by Hsieh, Bloch, and Larsen (1998). R² must lie in the interval [0, 1].

X distribution:

There a 7 options for the distribution of the X_i .

1. Binomial [P(k) = n!/(k! · (n - k)!) π^k (1 − π)^{N −k}, where k is the number of successes (X = 1) in N trials of a Bernoulli process with probability of success π, 0 < π < 1 ]
2. Exponential [f(x) = (1/λ)e^−1/λ , exponential distribution with parameter λ > 0]
3. Lognormal [f(x) = 1/(xσ√(2π)) exp[−(ln x − µ)²/(2σ²)], lognormal distribution with parameters µ and σ > 0.]
4. Normal [f(x)= 1/(σ√(2π)) exp[−(x − µ)²/(2σ²)], normal distribution with parameters µ and σ > 0)
5. Poisson (P(X = k) = (λ^k/k!)e^−λ, Poisson distribution with parameter λ > 0)
6. Uniform (f(x) = 1/(b − a) for a ≤ x ≤ b, f(x) = 0 otherwise, continuous uniform distribution in the interval [a, b], a < b)
7. Manual (allows to specify manually the variance of βˆ under H₀ and H₁)
   

G*Power provides two different types of procedure to calculate power: An enumeration procedure and large sample approximations. The Manual mode is available only in the large sample approximations.

Options

Input mode

You can choose between two input methods for the effect size: The effect size may be given by either specifying the two proportions p₁ and p₂ defined above, or instead by specifying p₁ and the odds ratio OR.

Procedure

G*Power provides two different types of procedure to estimate power. An enumeration procedure proposed by Lyles, Lin, and Williamson (2007) and large sample approximations.

The enumeration procedure seems to provide reasonably accurate results over a wide range of situations, but it can be rather slow and may need large amounts of computer memory. The large sample approximations are much faster.

Results of Monte-Carlo simulations indicate that the accuracy of the procedures proposed by Demidenko (2007) and Hsieh et al. (1998) are comparable to that of the enumeration procedure for N > 200. The procedure base on the work of Demidenko (2007) is more general and slightly more accurate than that proposed by Hsieh et al. (1998). We thus recommend to use the procedure proposed by Demidenko (2007) as the standard procedure.

The enumeration procedure of Lyles et al. (2007) may be used to validate the results (if the sample size is not too large). It must also be used if one wants to compute the power for likelihood ratio tests.

The enumeration procedure provides power analyses for the Wald test and the Likelihood ratio test. The general idea is to construct an exemplary data set with weights that represent response probabilities given the assumed values of the parameters of the X distribution. Then a fit procedure for the generalized linear model is used to estimate the variance of the regression weights (for Wald tests) or the likelihood ratio under H₀ and H₁ (for likelihood ratio tests).

The size of the exemplary data set increases with N and the enumeration procedure may thus be rather slow (and may need large amounts of computer memory) for large sample sizes. The procedure is especially slow for analysis types other then Post hoc power analyses, which internally call the power routine several times.
By specifying a threshold sample size N you can restrict the use of the enumeration procedure to sample sizes < N. For sample sizes ≥ N the large sample approximation selected in the Options dialog is used.

Note: If a computation takes too long you can abort it by pressing the ESC key.

G*Power provides two different large sample approximations for a Wald-type test. Both rely on the asymptotic normal distribution of the maximum likelihood estimator for parameter β₁ and are related to the method described by Whittemore (1981). The accuracy of these approximations increases with sample size, but the deviation from the true power may be quite noticeable for small and moderate sample sizes. This is especially true for X distributions that are not symmetric about the mean, i.e. the lognormal, exponential, and poisson distribution, and the binomial distribution with π ≠ 1/2.The approach of Hsieh et al. (1998) is restricted to binary covariates and covariates with a standard normal distribution. The approach based on Demidenko (2007) is more general and usually more accurate and is recommended as standard procedure. For this procedure, a variance correction option can be selected that compensates for variance distortions that may occur in skewed X distributions (see implementation notes). With the Hsieh procedure selected, the program automatically switches to the procedure of Demidenko if a distribution other than the standard normal or the binomial distribution is selected.

Examples

We first consider a model with a single continuous predictor X which is normally distributed with m = 0 and σ = 1. We assume that the event rate under H₀ is p₁ = 0.5 and that the event rate under H₁ is p₂ = 0.6. The odds ratio is then OR = (0.6/0.4)/(0.5/0.5) = 1.5, and we have β₁ = log(OR) ≈ 0.405. We want to estimate the sample size necessary to achieve, in a two-sided test with α = 0.05, a power of at least 0.95. We want to specify the effect size in terms of the odds ratio. Using the procedure of Hsieh et al. (1998) the input and output is as follows:

Select

Type of power analysis: A priori

Options:

Effect size input mode: Odds ratio
Procedure: Hsieh et al. (1998)

Input

Tail(s): Two
Odds ratio: 1.5
Pr(Y=1) H0: 0.5
α err prob: 0.05
Power (1-β err prob): 0.95
R2 other X: 0
X distribution: Normal
X parm µ: 0
X parm σ: 1

Output

Critical z: 1.959964
Total sample size: 317
Actual power: 0.950486

The results indicate that the necessary sample size is 317. This result replicates the value in Table II in Hsieh et al. (1998) for the same scenario. Using the other large sample approximation proposed by Demidenko (2007) with variance correction we arrive at N = 337. Without variance correction we arrive at N = 355.

In the enumeration procedure proposed by Lyles et al. (2007) the χ2 statistic is used and the output is

Noncentrality parameter λ: 13.029675
Critical χ2: 3.841459
Df: 1
Total sample size: 358
Actual power: 0.950498

Thus, this routine estimates the minimum sample size in this case to be N = 358.

In a Monte-Carlo simulation of the Wald test in the above scenario with 50000 independent cases we found a mean power of 0.940, 0.953, 0.962, and 0.963 for samples sizes 317, 337, 355, and 358, respectively. This indicates that in this case the method based on Demidenko (2007) with variance correction yields the best approximation.

We now assume that we have additional covariates and estimate the squared multiple population correlation with these others covariates to be R² = 0.1. All other conditions are identical. The only change we need to make is to set the input field R2 other X to 0.1. Using the procedure of Demidenko (2007) with variance correction under this condition we find that the necessary sample size increases from 337 to a value of 375.

As an example for a model with one binary predictor X we choose the values of the fourth example in Table I in Hsieh et al. (1998). That is, we assume that the event rate under H₀ is p₁ = 0.05, and the event rate under H₁ is p₂ = 0.1. We further assume a balanced design (π = 0.5) with equal sample frequencies for X = 0 and X = 1. Again we want to determine the sample size necessary to achieve in a two-sided test with α = 0.05 a power of at least 0.95. We want to specify the effect size directly in terms of p₁ and p₂:

Select

Type of power analysis: A priori

Options:

Effect size input mode: Two probabilities 
Procedure: Hsieh et al. (1998)

Input

Tail(s): Two
Pr(Y=1 |X=1) H1: 0.1
Pr(Y=1 |X=1) H0: 0.05
α err prob: 0.05
Power (1-β err prob): 0.95
R2 other X: 0
X Distribution: Binomial
X parm π: 0.5

Output

Critical z: 1.959964
Total sample size: 1437
Actual power: 0.950068

According to these results the necessary sample size is 1437. This replicates the sample size given in Table I in Hsieh et al. (1998) for this example.

This result is confirmed by the procedure proposed by Demidenko (2007) with variance correction. The procedure without variance correction and the procedure of Lyles et al. (2007) for the Wald test yield N = 1498. In Monte-Carlo simulations of the Wald test (50000 independent cases) we found mean power values of 0.953, and 0.961 for sample sizes of 1437 and 1498, respectively.

According to these results, the procedure of Demidenko (2007) with variance correction again yields the best power estimate for the tested scenario.

Let us select the procedure of Demidenko (2007) with variance correction in the Options dialog. Changing just the parameter of the binomial distribution (the prevalence rate) to a lower value of π = 0.2 increases the sample size to a value of 2158. Changing π to an equally unbalanced but higher value, 0.8, increases the required sample size further to 2368. These examples demonstrate the fact that a balanced design requires a smaller sample size than an unbalanced design, and a low prevalence rate requires a smaller sample size than a high prevalence rate (Hsieh et al., 1998, p. 1625).

Related tests

Poisson regression

Implementation notes

Enumeration procedure

The procedures for the Wald and Likelihood ratio tests are implemented exactly as described in Lyles et al. (2007).

Large sample approximations

The large sample procedures for the univariate case are both related to the approach outlined in Whittemore (1981). The correction for additional covariates has been proposed by Hsieh et al. (1998). As large sample approximations they get more accurate for larger sample sizes.

Demidenko procedure

In the procedure based on Demidenko (2007, 2008), the H₀ distribution is the standard normal distribution N(0, 1). The H₁ distribution the normal distribution N(m₁, s₁) with:

m₁ =(√(N(1 − R²)/v₁) β₁
s₁ = √((av₀ + (1 − a)v₁)/v₁)

where N denotes the sample size, R² the squared multiple correlation coefficient of the covariate of interest on the other covariates, and v₁ the variance of βˆ1 under H₁, whereas v₀ is the variance of ηˆ1 for H₀ with b^*₀ = ln(µ/(1 − µ)), with µ = ∫ fX(x) exp(b₀ + b₁ x)/(1 + exp(b₀ + b₁ x)dx. For the procedure without variance correction a = 0, that is, s₁ = 1. In the procedure with variance correction a = 0.75 for the lognormal distribution, a = 0.85 for the binomial distribution, and a = 1, that is, s₁ = √(v₀ /v₁) for all other distributions.

The motivation of the setting just described is as follows: Under H₀ the Wald-statistic has an asymptotic standard normal distribution, whereas under H₁, the asymptotic normal distribution of the test statistic has a mean of β₁/se(βˆ₁) = β₁/√(v₁/N), and a standard deviation of 1. However, the variance for finite n is biased (the degree depending on the X distribution), and √((av₀ + (1 − a)v₁)/v₁), which is close to 1 for symmetric X distributions but deviates from 1 for skewed distributions, often yields a much better estimate of the actual standard deviation of the distribution of βˆ₁ under H₁. This was confirmed in extensive Monte-Carlo simulations with a wide range of parameters and X distributions.

The procedure uses the result that the (m + 1) maximum likelihood estimators β₀, β₁, . . . , β_m are asymptotically normally distributed, where the variance-covariance matrix is given by the inverse of the (m + 1) × (m + 1) Fisher information matrix I. The (i, j)th element of I is given by



Thus, in the case of one continuous predictor, I is a 4 × 4 matrix with elements


where ƒX(x) is the PDF of the X distribution (for discrete predictors, the integrals must be replaced by corresponding sums). The element M₁₁ of the inverse of this matrix (M = I⁻¹), that is the variance of β₁, is given by: M₁₁ = Var(β) = I₀₀/(I₀₀ I₁₁ − I²₀₁). In G*Power, numerical integration is used to compute these integrals.

To estimate the variance of βˆ₁ under H₁, the parameters β₀ and β₁ in the equations for I_ij are chosen as implied by the input, that is, β₀ = log[p₁/(1 − p₁)], β₁ = log[OR]. To estimate the variance under H₀, one needs to choose β₁ = 0 and β₀ = β^*₀, where β^*₀ is chosen as defined above.

Hsieh et al. procedure

The procedures proposed in Hsieh et al. (1998) are used. The samples size formula for continuous, normally distributed predictors is [Eqn (1) in Hsieh et al. (1998)]:


where β* = log([p₁/(1 − p₁)]/[p₂/(1 − p₂)]) is the tested effect size, and p₁, p₂ are the event rates at the mean of X and one SD above the mean, respectively. For binary predictors the sample size formula is [see Eqn (2) in Hsieh et al. (1998)]:



where p = (1 − B)p₁ + B p₂ is the overall event rate, B is the proportion of the sample with X = 1, and p₁, p₂ are the event rates at X = 0 and X = 1 respectively.

Validation

To check the correct implementation of the procedure proposed by Hsieh et al. (1998), we replicated all examples presented in Tables I and II in Hsieh et al. (1998). The single deviation found for the first example in Table I on p. 1626 (sample size of 1281 instead of 1282) is probably due to rounding errors. Further checks were made against the corresponding routine in PASS (Hintze, 2006) and we usually found complete correspondence. For multiple logistic regression models with R2 other X > 0, however, our values deviated slightly from the result of PASS. We believe that our results are correct. There are some indications that the reason for these deviations is that PASS internally rounds or truncates sample sizes to integer values.

To validate the procedures of Demidenko (2007) and Lyles et al. (2007) we conducted Monte-Carlo simulations of Wald tests for a range of scenarios. In each case 150000 independent cases were used. This large number of cases is necessary to get about 3 digits precision. In our experience, the common praxis to use only 5000, 2000 or even 1000 independent cases in simulations (Lyles et al., 2007; Hsieh, 1989; Shieh, 2001) may lead to rather imprecise and thus misleading power estimates.

The table displayed below shows the errors in the power estimates for different procedures. Dem(c) and Dem denote the procedure of Demidenko (2007) with and without variance correction. LLW(W) and LLW(L) denote the procedure of Lyles et al. (2007) for the Wald test and the likelihood ratio test, respectively. All six predefined distributions were tested (the parameters are given in the table head). The following 6 combinations of Pr(Y=1|X=1) H0 and total sample size were used: (0.5,200),(0.2,300),(0.1,400),(0.05,600),(0.02,1000). These values were fully crossed with four odds ratios (1.3, 1.5, 1.7, 2.0), and two alpha values (0.01, 0.05). Max and mean errors were calculated for all power values < 0.999. The results show that the precision of the procedures depend on the X distribution. The procedure of Demidenko (2007) with the variance correction proposed here delivered the best predictions of the simulated power values.

References

Demidenko, E. (2007). Sample size determination for logistic regression revisited. Statistics in Medicine, 26, 3385-3397.

Demidenko, E. (2008). Sample size and optimal design for logistic regression with binary interaction. Statistics in Medicine, 27, 36-46.

Hintze, J. (2006). NCSS, PASS, and GESS. Kaysville, Utah: NCSS.

Hsieh, F. Y., Bloch, D. A., & Larsen, M. D. (1998). A simple method of sample size calculation for linear and logistic regression. Statistics in Medicine, 17, 1623-1634.

Lyles, R. H., Lin, H.-M., & Williamson, J. M. (2007). A practical approach to computing power for generalized linear models with nominal, count, or ordinal responses. Statistics in Medicine, 26, 1632-48.

Whittemore, A. S. (1981). Sample size for logistic regression with small response probability. Journal of the American Statistical Association, 76, 27-32.