SHARE
TWEET

Untitled

a guest Mar 17th, 2018 126 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # using factors (resulting in mean(Group)=0.5, mean(Sex)=0.5)
  2. sim = function(n=40, mu=0.6343) {
  3.   a  = rnorm(n,   0.0, 1)
  4.   b1 = rnorm(n/2, mu - mu/4, 1)
  5.   b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2
  6.  
  7.   dat = data.frame(Group = c(rep("C", n), rep("T", n)),
  8.                    Sex   = rep(c(rep("F", n/2), rep("M", n/2)), 2),
  9.                    Val   = c(a, b1, b2))
  10.   fit   = lm(Val ~ Group*Sex, dat)
  11.   return(summary(fit)$coefficients)
  12. }
  13.  
  14. # mean(Sex)=0, mean(Group)=0
  15. sim2 = function(n=40, mu=0.6343) {
  16.   a  = rnorm(n,   0.0, 1)
  17.   b1 = rnorm(n/2, mu - mu/4, 1)
  18.   b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2
  19.  
  20.   dat = data.frame(Group = c(rep(-0.5, n), rep(0.5, n)),
  21.                    Sex   = rep(c(rep(-0.5, n/2), rep(0.5, n/2)), 2),
  22.                    Val   = c(a, b1, b2))
  23.   fit   = lm(Val ~ Group*Sex, dat)
  24.   return(summary(fit)$coefficients)
  25. }
  26.  
  27. n = 40
  28. d = power.t.test(n=n, power=0.8)$delta
  29.  
  30. res=sapply(1:1000, function(x) sim(n, d)) # variant 1
  31. #res=sapply(1:1000, function(x) sim2(n, d)) # variant 2
  32.  
  33. mean(res[14,]<0.05) # real power main effect
  34. mean(res[16,]<0.05) # real power interaction
  35. data.frame(sd.of.coef.sample=sapply(1:4, function(x) sd(res[x,])),
  36.            mean.coef.se=sapply(5:8, function(x) mean(res[x,])),
  37.            q25.coef.se=sapply(5:8, function(x) quantile(res[x,], 0.25)),
  38.            q75.coef.se=sapply(5:8, function(x) quantile(res[x,], 0.75)),
  39.            coeff=sapply(1:4, function(x) mean(res[x,])),
  40.            z.coeff=sapply(1:4, function(x) mean(res[x,]/mean(res[x+4,]))))
  41.  
  42. mean(res[2, res[14,]<0.05]) # mean of significant main coefs
  43. mean(res[4, res[16,]<0.05]) # mean of significant interaction coefs
  44.  
  45. mean(res[2, res[14,]<0.05]) / d     # inflation main coef
  46. mean(res[4, res[16,]<0.05]) / (d/2) # inflation interaction coef
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top