-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathrsnmm.R
112 lines (111 loc) · 5.33 KB
/
rsnmm.R
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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
rsnmm <- function(n, tmax, control, ...) {
control <- if (missing(control)) rsnmm.control(...)
else do.call("rsnmm.control", control)
tmax <- tmax + (tmax %% 2) + 1
time <- rep(0:(tmax - 1), n)
tfun <- do.call("data.frame", lapply(control$tfun, function(f) f(time, tmax)))
coef.err <- 0
control$cormatrix <- matrix(control$coralpha, tmax, tmax)
diag(control$cormatrix) <- 1
if (control$corstr == "exchangeable") {
err <- sqrt(control$coralpha) * rep(rnorm(n, sd = control$sd), each = tmax)
err <- err + rnorm(n * tmax, sd = sqrt(with(control, sd^2 * (1 - coralpha))))
}
else {
## provisional error
err <- ifelse(time == 0, rnorm(n, sd = control$sd),
rnorm(n * (tmax - 1),
sd = sqrt(with(control, sd^2 * (1 - coralpha^2)))))
err[time == 0] <- rnorm(n, sd = control$sd)
coef.err <- control$coralpha
control$cormatrix <- matrix(with(control,
coralpha^(abs(row(cormatrix) -
col(cormatrix)))), tmax, tmax)
}
d <- .C("rsnmm",
n = as.integer(n),
tmax = as.integer(tmax),
ty = as.double(tfun$ty),
tmod = as.double(tfun$tmod),
tavail = as.double(tfun$tavail),
tstate = as.double(tfun$tstate),
beta = with(control, as.double(c(beta0, beta1))),
eta = as.double(control$eta),
mu = as.double(control$mu),
theta = with(control, as.double(c(theta0, theta1))),
coef.avail = as.double(control$coef.avail),
coef.state = as.double(control$coef.state),
coef.err = as.double(coef.err),
avail = as.integer(rep(0, n * tmax)),
base = as.double(rep(rnorm(n), each = tmax)),
state = as.integer(rep(0, n * tmax)),
a = as.integer(rep(0, n * tmax)),
prob = as.double(rep(0, n * tmax)),
y = as.double(rep(0, n * tmax)),
err = as.double(err),
state.center = as.double(rep(0, n*tmax)),
a.center = as.double(rep(0, n*tmax)),
avail.center = as.double(rep(0, n*tmax)))
d <- data.frame(id = rep(1:n, each = tmax), time = time,
ty = d$ty, tmod = d$tmod, tavail = d$tavail, tstate = d$tstate,
base = d$base, state = d$state, a = d$a, y = d$y, err = d$err,
avail = d$avail, prob = d$p, a.center = d$a.center,
state.center = d$state.center, avail.center = d$avail.center,
one = 1)
## nb: for a given row, y is the proximal response
d$lag1y <- with(d, delay(id, time, y))
d$lag2y <- with(d, delay(id, time, y, 2))
d$lag1err <- with(d, delay(id, time, err))
d$lag1avail <- with(d, delay(id, time, avail))
d$lag1avail.center <- with(d, delay(id, time, avail.center))
d$lag2avail <- with(d, delay(id, time, avail, 2))
d$lag2avail.center <- with(d, delay(id, time, avail.center, 2))
d$lag1a <- with(d, delay(id, time, a))
d$lag2a <- with(d, delay(id, time, a, 2))
d$lag1prob <- with(d, delay(id, time, prob))
d$lag2prob <- with(d, delay(id, time, prob, 2))
d$lag1a.center <- with(d, delay(id, time, a.center))
d$lag2a.center <- with(d, delay(id, time, a.center, 2))
d$lag1tmod <- with(d, delay(id, time, tmod))
d$lag2tmod <- with(d, delay(id, time, tmod, 2))
d$lag1state <- with(d, delay(id, time, state))
d$lag1state.center <- with(d, delay(id, time, state.center))
rownames(d) <- NULL
attributes(d) <- c(attributes(d), control)
d
}
rsnmm.control <- function(origin = 1, sd = 1,
coralpha = sqrt(0.5),
corstr = c("ar1", "exchangeable"),
beta0 = c(-0.2, 0, 0, 0.2, 0), beta1 = rep(0, 4),
eta = c(0, 0, 0.8, -0.8, 0), mu = rep(0, 3),
theta0 = c(0, 0.8), theta1 = c(0, 0),
coef.avail = c(100, rep(0, 3)), coef.state = rep(0, 5),
tfun = NULL, lag = 3 + any(beta1 != 0)) {
corstr <- match.arg(corstr)
if (is.null(tfun))
tfun <- rep(list(function(tcur, tmax) rep(0, length(tcur))), 4)
list(origin = 1, lag = lag,
## error SD, correlation
sd = sd, coralpha = coralpha, corstr = corstr,
## proximal effect coefficients
beta0 = setNames(beta0, c("one", "tmod", "base", "state", "lag1a")),
## delayed effect coefficients
beta1 = setNames(beta1, c("one", "lag1tmod", "base", "lag1state")),
## treatment probability model coefficients
eta = setNames(eta, c("one", "base", "state", "lag1a", "lag1y")),
## exogenous or time-invariant main effects
mu = setNames(mu, c("one", "ty", "base")),
## time-varying main effects, centered and proximal
theta0 = setNames(theta0, c("avail", "state")),
## time-varying main effects, centered and delayed
theta1 = setNames(theta1, c("lag1avail", "lag1state")),
## availability model coefficients
coef.avail = setNames(coef.avail, c("one", "tavail", "lag1a", "lag1y")),
## binary state model coefficients
coef.state = setNames(coef.state,
c("one", "tstate", "base", "lag1state", "lag1a")),
## functions of time in the main effect, proximal effect,
## availability model, and binary state model
tfun = setNames(tfun, c("ty", "tmod", "tavail", "tstate")))
}