1

Although there has been some detailed discussions about power analysis on this website (for example here and here), the answer provided to this question has outlines the steps to simulating a power analysis, here.

Say we take some data (data was linked to a bootstrapping question)

We create a regression that will predict admit based on the two continous variables gpa and gre

  • Now we have a n=400.
  • We can then elect our power level, alpha = 0.5
  • The effect size you would like to detect, e.g., odds ratios (we obtain this from our regression)

So in following the detailed method provided by @gung here, I want to run the simulation. Here is the code I have adjusted, but my output is not correct. Can someone outline what I have not understood

mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
head(mydata)

set.seed(1234)

my.mod <- glm(admit ~ gre + gpa , data = mydata, family = "binomial")


repetitions <- length(mydata$admit)

gre <- mydata$gre
gpa <- mydata$gpa


significant = matrix(nrow=repetitions, ncol=4)

for(i in 1:repetitions){
  responses          = mydata$admit
  #responses          = rbinom(n=N, size=1, prob=mydata$admit)      # we can interchange this comment
  model              = glm(responses ~ gre + gpa, family = binomial(link="logit"))
  significant[i,1:2] = (summary(model)$coefficients[2:3,4]<.05)
  significant[i,3]   = sum(significant[i,1:2])
  modelDev           = model$null.deviance-model$deviance
  significant[i,4]   = (1-pchisq(modelDev, 2))<.05
}



sum(significant[,1])/repetitions      # pre-specified effect power for gre

sum(significant[,2])/repetitions      # pre-specified effect power for gpa

sum(significant[,4])/repetitions  # power for likelihood ratio test of model

sum(significant[,3]==2)/repetitions   # all effects power

sum(significant[,3]>0)/repetitions    # any effect power
Community
  • 1
user08041991
  • 255
  • 3
  • 6
  • It looks like you haven't actually simulated anything; you've just done the same thing multiple times to the same data. Even if you use `responses = rbinom(n=N, size=1, prob=mydata$admit)`, that's not correct because `admit` is the dependent variable (taking values of 0 or 1) rather than being probabilities. – mark999 Feb 27 '16 at 20:49

0 Answers0