Guest User

Untitled

a guest
Mar 17th, 2018
132
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. sim1 = function(n=40, mu=0.6343) {
  2.   a  = rnorm(n,   0.0, 1)
  3.   b1 = rnorm(n/2, mu - mu/4, 1)
  4.   b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2
  5.  
  6.   dat = data.frame(Group = c(rep("C", n), rep("T", n)),
  7.                    Sex   = rep(c(rep("F", n/2), rep("M", n/2)), 2),
  8.                    Val   = c(a, b1, b2))
  9.   fit   = lm(Val ~ Group*Sex, dat)
  10.   return(summary(fit)$coefficients)
  11. }
  12.  
  13. sim2 = function(n=40, mu=0.6343) {
  14.   a  = rnorm(n,   0.0, 1)
  15.   b1 = rnorm(n/2, mu - mu/4, 1)
  16.   b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2
  17.  
  18.   dat = data.frame(Group = c(rep("C", n), rep("T", n)),
  19.                    Sex   = rep(c(rep("M", n/2), rep("F", n/2)), 2),  # we only exchange labels M & F
  20.                    Val   = c(a, b1, b2))
  21.   fit   = lm(Val ~ Group*Sex, dat)
  22.   return(summary(fit)$coefficients)
  23. }
  24.  
  25. n = 40
  26. d = power.t.test(n=n, power=0.8)$delta
  27.  
  28. # variation 1
  29. res1=sapply(1:1000, function(x) sim1(n, d))
  30. mean(res1[14,]<0.05) # power main effect
  31.  
  32. # variation 2
  33. res2=sapply(1:1000, function(x) sim2(n, d))
  34. mean(res2[14,]<0.05) # power main effect
RAW Paste Data