# Power analysis for probit model in R

New to power analysis,

I am trying to perform a power analysis in R to determine n individuals needed to achieve 0.80 power for a probit regression.

I've tried looking into prw but there doesn't seem to have any functions for GLMs. The closest option that I see is pwr.f2.test: Power calculations for the general linear model But this seems inappropriate because my data is not normally distributed.

I've found some other packages such as SIMR that are used for generalized linear mixed models.

Edit: The probit model includes a binary outcome and a single continuous ordinal predictor. The predictor is different exposure times of a solution. I’m trying to determine the lethal time of exposure (LT50 etc). I’m trying to determine how many individuals I would need for the model to fit the data well (probably a gof?).

My question is: What is the right way to perform a power analysis in this scenario?

• What does your probit model look like? What else is in the model besides the binary outcome? Jan 13, 2023 at 2:25
• @dimitriy My probit model has a single continuous (ordinal) predictor and with the binary outcome. I’m trying to determine Lethal time of a treatment. So the continuous predictor is different exposure times of the same concentration of a solution. So, I’m looking to determine how many individuals I would need for the model to for the data well (gof test, I believe). Jan 13, 2023 at 12:00

I have not seen a probit sample size calculation like this, but I have seen this done for a logit. I suspect the choice of model will not matter substantially unless you are way out in the tails (very rare or very frequent events). This could be verified with a simulation.

You will need to know or guesstimate the probability of success at the mean of $$x$$ and at one standard deviation above the mean of $$x$$: \begin{align} p_1 &= \Pr(y=1 \vert x= \bar x) \\ p_2 &= \Pr(y=1 \vert x= \bar x + \sigma_x)\end{align} You can usually get these from a pilot study or from previous work.

For example,

> ss <- powerlog(p1 = 0.5, p2 = 0.6, alpha = .05, power = 0.8, help = TRUE)

Logistic regression power analysis
One-tailed test: alpha = 0.05 , power = 0.8 , p1 = 0.5 , p2 = 0.6 , r2 = 0 , odds ratio = 1.5

n
164

Explanation of terms:
p1 is the probability that the outcome equals 1 when the predictor x is at the mean
p2 is the probability that the outcome equals 1 when the predictor x is one standard deviation above its mean
n is the sample size required for the logit model


Here is an R function that produced the output above:

# Logistic Regression Sample Size Analysis
#
# This function computes the sample size required for a logistic regression with a single continuous covariate of interest x
# and one or more optional control variables  z1-zk, Pr(y=1 | x, z1,...zk),
#
# It is based on Phil Ender's powerlog.ado (version 1.0) for Stata, which is itself based on POWERLOG.SAS (version 1.1)
# for SAS by Michael Friendly.
#
# Inputs:
# p1 is the probability that the outcome equals 1 when the predictor x is at its mean.
# p2 is the probability that the outcome equals 1 when the predictor x is one standard deviation above its mean.
# Both p1 and p2 must be strictly greater than 0 and strictly less than 1. The best way to obtain these numbers is through a
# pilot study that closely mimics your research design or from previous literature.
# r2 is the squared multiple correlation between the predictor variable x and all other variables in the model z1-zk (default = 0).
# Use r2=0 for a one-predictor model with only x.
# alpha is the significance level for the test. It is the probability of rejecting the null hypothesis when the null hypothesis is true (default = 5%)
# power is the probability of detecting a “true” effect when it exists, defined as 1 minus the probability of incorrectly accepting the null hypothesis (default = 80%)
# If help is not set to FALSE, the code returns an explanation (default = TRUE)
#
# Output:
# Prints the parameters of the hypothesis test and required sample size n
# If help is set to TRUE, prints out an explanation
# Returns a list containing p1, p2, r2, alpha, power, and n
#
# Example:
# ss <- powerlog(p1 = 0.25, p2 = 0.27, alpha = .05, power = 0.8, help = TRUE)
#
# Example Output:
# Logistic regression power analysis
# One-tailed test: alpha = 0.05 , power = 0.8 , p1 = 0.25 , p2 = 0.27 , r2 = 0 , odds ratio = 1.109589
# n
# 3,440
#
# Explanation of terms:
# p1 is the probability that the outcome equals 1 when the predictor x is at the mean
# p2 is the probability that the outcome equals 1 when the predictor x is one standard deviation above its mean
# r2 is the squared multiple correlation between the predictor variable and all other variables in the model
# n is the sample size required for the logit model
#

powerlog <- function(p1, p2, r2=0, alpha=0.05, power = 0.8, help= TRUE) {
pd <- p2 - p1
l1 <- p1 / (1-p1)
l2 <- p2 / (1-p2)
theta <- l2 / l1
or <- theta
lambda <- log(theta)
lambda2 <- lambda^2
za <- qnorm(1-alpha)

cat("\nLogistic regression power analysis\n")
cat("One-tailed test: alpha =", alpha, ", power =", power, ", p1 =", p1, ", p2 =", p2, ", r2 =", r2, ", odds ratio =", or, "\n\n")
cat("n\n")
delta <- (1 + (1 + lambda2) * exp(5 * lambda2 / 4)) / (1 + exp(-1 * lambda2 / 4))
zb <- qnorm(power)
N <- ((za + zb * exp(-1 * lambda2 / 4))^2 * (1 + 2 * p1 * delta)) / (p1 * lambda2)
N <- round(N / (1-r2))
cat(formatC(N, format="d", big.mark=","), "\n")

if (help != FALSE) {
cat("\nExplanation of terms:\n")
cat("p1 is the probability that the outcome equals 1 when the predictor x is at the mean\n")
cat("p2 is the probability that the outcome equals 1 when the predictor x is one standard deviation above its mean\n")
if (r2 != 0) {
cat("r2 is the squared mulitple correlation between the predictor variable and all other variables in the model\n")
}
cat("n is the sample size required for the logit model\n")
}

return(list(p1 = p1, p2 = p2, r2 = r2, alpha = alpha, power = power, n = N))

}


For completeness, there is also an R package called WebPower that can handle logistic SS calculations. It also has a web interface. It uses a different approach and assumptions than I do, and it can do both two-sided and one-sided tests. It gets 169 rather than 164:

> WebPower::wp.logistic(n = NULL, p0 = 0.5, p1 = 0.6, alpha = 0.05,
+             power = 0.8, family = "normal", parameter = c(0,1), alternative = "greater")
Power for logistic regression

p0  p1 beta0     beta1       n alpha power
0.5 0.6     0 0.4054651 168.776  0.05   0.8


Response to Question in Comments

Here's a simplified example of how you might do a calculation like this. Say you are a ride-hailing company like Uber, and you want to estimate the probability of ordering a ride ($$y=1$$) as a function of price $$x$$ using an experiment where you vary the price that riders see in the app.

But say these are the early days, and you don't have a model for $$\Pr(\mathtt{ride} \vert x)$$. After all, the point of the experiment is to estimate this model reliably with the experiment you are designing. What is to be done? You might know something about the $$\Pr(\mathtt{ride} \vert x_1)$$, where $$x_1$$ is the current price. Taking the average over all the potential current customers that open your app and see a price would give you a $$\hat p_1$$. If you don't have a current price, you can use a focus group to estimate $$p_1$$ for a price in the middle of the range you are considering.

You should also have some idea of the range of $$x$$ that you want to test. Let's assume that range is from $$0$$ (free) to $$2 \cdot x_1$$ (2X current price), where each price is equally likely. So $$x$$ has a uniform distribution over $$[0,2 \cdot x_1]$$. The standard deviation of a uniform RV like that is $$\sigma_x = \sqrt(\frac{4}{12}) \cdot x_1$$. You do a pilot with price set to $$x_2 = x_1 + \sqrt(\frac{4}{12}) \cdot x_1$$ and get the average probability of sale $$\hat p_2$$.

If you can't do a pilot, you can make an assumption like this. In previous research, the short-run price elasticity of taxis was found to be $$0.6$$, so for a $$1\%$$ increase in price, quantity demanded falls by $$0.6\%$$. An increase from $$x_1$$ to $$x_2$$ is about $$60\%$$, so you would expect that $$\hat p_2 = 0.35 \cdot \hat p_1$$.

This gives you $$p_1$$ and $$p_2$$ to plug into the sample size calculator above.

• Hi @dimitriy, thank you so much for this! I do have a question regarding determining the probability of success at the mean of x. Would I simply plug the mean into the the model's equation to get the probability of an event occurring at the mean? Jan 19, 2023 at 21:33
• @scott.pilgrim.vs.r I added a toy example above since it would not fit into the comments. Jan 20, 2023 at 0:22
• That helps a lot! Thanks you! Jan 20, 2023 at 1:40