forked from schoam4/mc_power_med
-
Notifications
You must be signed in to change notification settings - Fork 0
/
.Rhistory
88 lines (88 loc) · 2.56 KB
/
.Rhistory
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
shiny::runApp()
runApp()
corMat <- diag(4)
corMat[1,2] <- corMar[2,1] <- .5
corMat[1,2] <- corMat[2,1] <- .5
corMat
corMat[2,4] <- corMat[4,2] <- .5
corMat
covMat <- corMat
powRep <- function(seed = 1234, Ns = N, covMatp = corMat){
set.seed(seed)
require(MASS)
incProgress(1/powReps)
dat <- mvrnorm(Ns, mu = c(0,0,0,0), Sigma = covMatp)
# Run regressions
m1 <- lm(dat[,2] ~ dat[,1])
m2 <- lm(dat[,3] ~ dat[,1] + dat[,2])
m3 <- lm(dat[,4] ~ dat[,2] + dat[,3] + dat[,1])
# Output parameter estimates and standard errors
a1 <- rnorm(mcmcReps, coef(m1)[2], sqrt(vcov(m1)[2,2]))
a2 <- rnorm(mcmcReps, coef(m2)[2], sqrt(vcov(m2)[2,2]))
b1 <- rnorm(mcmcReps, coef(m3)[2], sqrt(vcov(m3)[2,2]))
b2 <- rnorm(mcmcReps, coef(m3)[3], sqrt(vcov(m3)[3,3]))
d <- rnorm(mcmcReps, coef(m2)[3], sqrt(vcov(m2)[3,3]))
a1b1 <- a1*b1
a2b2 <- a2*b2
a1db2 <- a1*d*b2
# Calculate confidence intervals
low <- (1 - (conf / 100)) / 2
upp <- ((1 - conf / 100) / 2) + (conf / 100)
LL1 <- quantile(a1b1, low)
UL1 <- quantile(a1b1, upp)
LL2 <- quantile(a2b2, low)
UL2 <- quantile(a2b2, upp)
LLd <- quantile(a1db2, low)
ULd <- quantile(a1db2, upp)
# Is rep significant?
c(LL1*UL1 > 0, LL2*UL2 > 0, LLd*ULd > 0)
}
N = 200
powRep()
powRep <- function(seed = 1234, Ns = N, covMatp = corMat){
set.seed(seed)
require(MASS)
#incProgress(1/powReps)
dat <- mvrnorm(Ns, mu = c(0,0,0,0), Sigma = covMatp)
# Run regressions
m1 <- lm(dat[,2] ~ dat[,1])
m2 <- lm(dat[,3] ~ dat[,1] + dat[,2])
m3 <- lm(dat[,4] ~ dat[,2] + dat[,3] + dat[,1])
# Output parameter estimates and standard errors
a1 <- rnorm(mcmcReps, coef(m1)[2], sqrt(vcov(m1)[2,2]))
a2 <- rnorm(mcmcReps, coef(m2)[2], sqrt(vcov(m2)[2,2]))
b1 <- rnorm(mcmcReps, coef(m3)[2], sqrt(vcov(m3)[2,2]))
b2 <- rnorm(mcmcReps, coef(m3)[3], sqrt(vcov(m3)[3,3]))
d <- rnorm(mcmcReps, coef(m2)[3], sqrt(vcov(m2)[3,3]))
a1b1 <- a1*b1
a2b2 <- a2*b2
a1db2 <- a1*d*b2
# Calculate confidence intervals
low <- (1 - (conf / 100)) / 2
upp <- ((1 - conf / 100) / 2) + (conf / 100)
LL1 <- quantile(a1b1, low)
UL1 <- quantile(a1b1, upp)
LL2 <- quantile(a2b2, low)
UL2 <- quantile(a2b2, upp)
LLd <- quantile(a1db2, low)
ULd <- quantile(a1db2, upp)
# Is rep significant?
c(LL1*UL1 > 0, LL2*UL2 > 0, LLd*ULd > 0)
}
powRep()
mcmcReps = 20000
powRep()
conf = 95
powRep()
pow <- lapply(sample(1:50000, powReps), powRep)
powReps = 500
pow <- lapply(sample(1:50000, powReps), powRep)
colSums(matrix(unlist(pow), nrow = powReps)) / powReps
unlist(pow)
df <- data.frame("Parameter" = c("a1b1", "a2b2", "a1db2"),
"N" = rep(N, 3),
"Power" = colSums(matrix(unlist(pow), nrow = powReps, byrow = TRUE)) / powReps)
df
runApp()
runApp()
runApp()