# 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