R - Power analysis for a logistic regression

Question

I am looking for a way to make a proper power analysis for a logistical regression in R. So, I created this model for a logistical regression:

model <- glm(df1$fish. ~ df1$Temperature, family="binomial", data=df1)

options(scipen=999)
summary(model)

# McFadden’s R2 to test if model fits
library("pscl")
pscl::pR2(model)["McFadden"]

To evaluate if I have a good sample size, I wanted to make a power analysis. Here

https://cran.r-project.org/web/packages/pwr/vignettes/pwr-vignette.html?fbclid=IwAR3VJR1zN2cd0UgE9xzKOQRgjthv-UFQoxXUBrZa8G8FfnWfBHbeR-9kLU8

I found a code to use for that.

library("pwr")
plot(pwr.f2.test(u=1, v=(1:48), f2=?, sig.level=0.05)$power,
     3:50, xlab="power", ylab="sample size") 
# u = number of coefficients in the model minus the intercept (2 - 1 = 1)
# v = number of error degrees of freedom (v = n - u - 1)
# f2 = effect size = ?
abline(h=19) # enter my sample size 
pwr.f2.test(u=1, v=(19-1-1), f2=?, sig.level=0.05)$power # output number (percentage probability)
However, my problem now is, that I don't know what to enter for an effect size f2. On the link above it says, that you use R2/(1-R2). In my case, I don't have an R2, only a MC Fadden's R2 (~0.495555), which is obviously not the same, so I cannot use it the same way.

Can someone explain to me, what exactly the McFadden's R2, as well as the f2 is and what it tells me? Does anyone know what I have to enter for f2 in my case?

I used the following steps, to validate my regression model:

# make predictions
predicted <- predict(model, df1, type="response")
predicted
# predicts the probability of fish appearing at a certain temperature

# validate model with confusion matrix
confmatrix <- table(actualValue = df1$fish, predictedValue = predicted > 0.5)
confmatrix
# percentage accuracy
(confmatrix[[1, 1]] + confmatrix[[2, 2]]) / sum(confmatrix)

#find optimal cutoff probability to use to maximize accuracy
library(InformationValue)
optimal <- optimalCutoff(df1$fish, predicted)[1]
optimal

# calculate sensitivity: “true positive rate”
sensitivity(df1$fish, predicted)

#calculate specificity: “true negative rate”
specificity(df1$fish, predicted)

# calculate total misclassification error rate
misClassError(df1$fish, predicted, threshold=optimal)

#plot the ROC curve
plotROC(df1$fish, predicted)

My thought was, whether I can use the confmatrix as an effect size? It seems to be similar to R2, however I am not sure if I have understood it correctly. It is just an idea I had, but I am only guessing here.

Answer 1

Edited Post

Before going any further, I should mention that you seem still unclear of what you want to do. Indeed, you have several questions in your post mobilizing very different concepts. I will try to address them as clear as possible, presenting different use case:

1. Just to clarify one thing, power analysis refers to as simulation studies. In this case you want to assess whether your model is good to detect signal (rejecting the null hypothesis). Typically, we will simulate a large number of replicates, playing on different parameters. Traditionally, we like to play with beta parameters (the association level between variables and my outcome), number of variables, nature of variables, and number of individuals. It is worth noting that the number of replicates is different from your number of individuals. The former is independent of your model (it permits you to test the "asymptotic" performance of your model) while the latter directly depends on your model. In your context, "How many individuals I need to have a certain level of power?" for example. Here your power level is not specified a priori. However, when you want to have the "optimal" sample size, you must specify your target power. In practice we want power greater than 0.8 but it depends on your context.

2. If I am correct you cannot use pwr package. It was developed only for linear models. Others packages can do power analyses for logistic regressions. Please be aware of the hypotheses (continuous predictors for example).

3. Finally, here it depends on what you want to do. I provided a brief example to illustrate how to do power analysis with logistic regression exploiting the different notions you mentioned in your post. In my code I illustrated very simple applications of simulation studies, such as, power (how many I rejected the null association), bias calculation (How my estimates differs from my ground truth (true values)?) or model performance. For the latter you can use any metric you believe is relevant in your context.

#Let's simulate data by ourselves
#We simulated 1000 replicates and we want to assess power (Beta1 < alpha), bias, or computing 
#any statistic

#Assuming these values for the log odds ratio
beta0= 0.5
beta1= 0.75

#Your significance threshold
alpha=0.05

#Number of replicates 
n = 1000

#If we are interested in rejecting the null
power.vec = c()

#If we are interested in assessing parameter bias
bias.vec = c()

#If we are interested in computing Mac Fadden pseudo R2 
R2.vec = c()

for(i in 1:n){
  
  #Assuming one continous predictor and 100 individuals
  X = rnorm(100)
  
 
  #Our logistic model
  pr = exp(beta0+X*beta1)/(1+exp(beta0+X*beta1))
  
  y = rbinom(100, 1,pr)
  
  model = glm(y~X, family = binomial(link="logit"))
  sm = summary(model)
  
  power.vec = c(power.vec, sm$coefficients[2,4]<alpha)
  bias.vec = c(bias.vec, beta1-sm$coefficients[2,1])
  R2.vec = c(R2.vec, 1-(sm$deviance/sm$null.deviance))
}

#The power (here rejecting the null for beta 1 depends on the number of individuals)
sum(power.vec)/n

#The bias is reducing when increasing the number of individuals
mean(abs(bias.vec))

#Mean pseudo R2
mean(R2.vec)

#Sample size calculation 

library(WebPower)

#p0 is probability of Y = 1 when X = 0: You should change this parameter
#p1 is probability of Y=1 when X=1: You should change this parameter

#alpha: your significance threshold 
#power is the target power
#family: the nature of your predictor

WebPower::wp.logistic(n=NULL, p0=0.10, p1=0.20, alpha=0.05, power=0.8, alternative = "two.sided", family="normal")

To conclude you can play on every parameters, e.g., number of individuals, effect measure, number of predictors, type of predictors etc....

For sample size calculations, I added an example from: https://med.und.edu/research/daccota/_files/pdfs/berdc_resource_pdfs/sample_size_r_module.pdf. Reducing the target power for the same probabilities and the same significance threshold will lead to smaller needed sample sizes. I will let you play with this.

Hope this helps.