0
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ So I got the hint, that for the logical regression, I should probably rather use the code that is suggested here: rdrr.io/cran/WebPower/man/wp.logistic.html. However, then I am still quite inexperienced with R and don't know exactly what to enter for p0 and p1. Is it the same as sensitivity(“true positive rate”) and specificity(“true negative rate”) and can be used analogue? And what should I use for 'alternative' and 'family'? Can I leave all suggestions that are given in the code, or do I have to chose, i.e. 'Bernoulli'? If I have to chose, which one would be necessary in my case? $\endgroup$
    – Feli
    Feb 9, 2023 at 21:06
  • $\begingroup$ Example and code here. $\endgroup$
    – dimitriy
    Feb 9, 2023 at 22:11
  • $\begingroup$ @dimitriy: I find this link helpful, however, do you think you are able to tell me how to calculate or figure out p0 and p1? I tried some things that didn't work and am past my limit of experience now $\endgroup$
    – Feli
    Feb 10, 2023 at 13:02
  • $\begingroup$ I am not sure if I can say anything else beyond what I have written. I would start another question, where you show your code, its results and explain what exactly did not work. If it's programming, that is probably best for stackoverflow. If it’s statistics, CV is the right place. $\endgroup$
    – dimitriy
    Feb 10, 2023 at 17:32

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Thank you so much for the detailed answer. However, since I am only starting to crawl with R and statistics in general, I am having difficulties understanding all of the variables you used. Can you elaborate on what exactly beta0 and beta1 signify and how I chose the values to enter there? Furthermore, what is the number of replicas as opposed to the number of individuals? What does the bias tell me about my data and is there a way to apply the power analysis to my logistic regression code I have above? I apologize, but I struggle to translate your approach to my data. $\endgroup$
    – Feli
    Feb 10, 2023 at 12:58
  • $\begingroup$ @Feli I edited the post to match your questions. $\endgroup$ Feb 10, 2023 at 15:38
  • $\begingroup$ Thank you so much, this was a big help, I think I finally found what I was looking for! $\endgroup$
    – Feli
    Feb 10, 2023 at 16:28
  • $\begingroup$ @Feli My pleasure! Please accept the answer if I answered your initial post and following questions. $\endgroup$ Feb 10, 2023 at 16:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.