• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Mar 17th, 2018 143 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.
Not a member of Pastebin yet?