Guest User

Untitled

a guest
Mar 17th, 2018
163
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