Nothing
context('multivar')
## Generate test data without littering the environment with temporary
## variables
x <- NULL
y <- NULL
d <- NULL
o <- NULL
local({
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
N <- 3
T <- 2
K <- 2
beta <- c(.1, -.5)
dd <- generate_data(N=N, T=T, K=K, beta=beta)
x <<- dd$x
y <<- dd$y
d <<- data.frame(i = rep(seq_len(N), each=T+1),
t = rep(seq_len(T+1), N),
matrix(aperm(x, c(1, 3, 2)), N*(T+1), K,
dimnames = list(NULL, paste0('x', seq(K)))),
y = c(y))
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
o <<- opm(y~x1+x2, d, n.samp = 10)
})
## Sanity check
test_that('data', {
expect_equal(x, array(c(-1.31047564655221, -0.289083794010798,
-0.349228549405948, -0.679491608575424,
-1.19566197009996, 1.03691313680308,
-0.23017748948328, -1.26506123460653,
0.11068271594512, 0.129287735160946,
1.22408179743946, 0.497850478229239,
2.30870831414912, 0.0631471481064739,
0.194158865245925, 2.46506498688328,
1.10981382705736, -1.21661715662964),
dim = c(3, 2, 3)))
expect_equal(y, matrix(c(-0.0899458588038237, -0.694025238411308,
-2.52543131039704, -0.560453024256735,
-2.04477798261599, -2.94693926957052,
-1.06948536801357, -0.812226112015961,
2.05939845332596),
nrow = 3, ncol = 3))
})
test_that('default', {
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
expect_equal(opm(x, y, n = 10),
structure(list(samples = list(rho = c(0.805, 0.425,
-0.181, 0.763,
0.581, 0.965,
0.0559999999999999, 0.533,
0.648, -0.086),
sig2 = 1 / c(11.6175662271191,
0.919115793520271,
0.989369792715488,
0.00440918521157215,
1.35621264272223,
4.24440799476086,
0.0931249250470442,
1.56847551750474,
0.371922970545138,
5.00022577787265),
beta = matrix(c(0.149966256134789, -1.07201170366086,
-3.32784640301039, 6.31976053552415,
0.399370400131068, -0.409833517459282,
0.942409502841064, -1.01613441888489,
-1.06148590924992, -1.70364397749032,
-0.746392848950218, -1.22143920569537,
-1.63866906134468, 7.01947247537598,
-0.536929048279301, -0.845295844019301,
-0.565341215754439, 0.0765501718682319,
-0.171463431207414, -1.13494915857112),
10, 2)),
fitted.values = matrix(c(0.878931188971585,
-0.878931188971585, 0.665670601363721, -0.665670601363721, -0.787856768435423,
0.787856768435423), 2, 3),
residuals = matrix(c(0.0367718470212829,
-0.0367718470212828, -0.214589957886459, 0.214589957886458, -0.647955514235535,
0.647955514235536), 2, 3),
df.residual = 2*3 - 3,
logLik = 0.704787094195681,
design = "balanced",
call = quote(opm(x = x, y = y, n.samp = 10)),
.Environment = environment(),
time.indicators = FALSE),
class = 'opm'))
})
test_that('named matrix', {
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
colnames(x) <- c('x1', 'x2')
expect_equal(opm(x, y, n = 10),
structure(list(samples = list(rho = c(0.805, 0.425,
-0.181, 0.763,
0.581, 0.965,
0.0559999999999999, 0.533,
0.648, -0.086),
sig2 = 1 / c(11.6175662271191,
0.919115793520271,
0.989369792715488,
0.00440918521157215,
1.35621264272223,
4.24440799476086,
0.0931249250470442,
1.56847551750474,
0.371922970545138,
5.00022577787265),
beta = matrix(c(0.149966256134789, -1.07201170366086,
-3.32784640301039, 6.31976053552415,
0.399370400131068, -0.409833517459282,
0.942409502841064, -1.01613441888489,
-1.06148590924992, -1.70364397749032,
-0.746392848950218, -1.22143920569537,
-1.63866906134468, 7.01947247537598,
-0.536929048279301, -0.845295844019301,
-0.565341215754439, 0.0765501718682319,
-0.171463431207414, -1.13494915857112),
10, 2, dimnames = list(NULL, c('x1', 'x2')))),
fitted.values = matrix(c(0.878931188971585,
-0.878931188971585, 0.665670601363721, -0.665670601363721, -0.787856768435423,
0.787856768435423), 2, 3),
residuals = matrix(c(0.0367718470212829,
-0.0367718470212828, -0.214589957886459, 0.214589957886458, -0.647955514235535,
0.647955514235536), 2, 3),
df.residual = 2*3 - 3,
logLik = 0.704787094195681,
design = "balanced",
call = quote(opm(x = x, y = y, n.samp = 10)),
.Environment = environment(),
time.indicators = FALSE),
class = 'opm'))
})
test_that('formula', {
expected <- structure(list(samples = list(rho = c(0.805, 0.425,
-0.181, 0.763,
0.581, 0.965,
0.0559999999999999, 0.533,
0.648, -0.086),
sig2 = 1 / c(11.6175662271191,
0.919115793520271,
0.989369792715488,
0.00440918521157215,
1.35621264272223,
4.24440799476086,
0.0931249250470442,
1.56847551750474,
0.371922970545138,
5.00022577787265),
beta = matrix(c(0.149966256134789, -1.07201170366086,
-3.32784640301039, 6.31976053552415,
0.399370400131068, -0.409833517459282,
0.942409502841064, -1.01613441888489,
-1.06148590924992, -1.70364397749032,
-0.746392848950218, -1.22143920569537,
-1.63866906134468, 7.01947247537598,
-0.536929048279301, -0.845295844019301,
-0.565341215754439, 0.0765501718682319,
-0.171463431207414, -1.13494915857112),
10, 2, dimnames = list(NULL, c('x1', 'x2')))),
fitted.values = matrix(c(0.878931188971585,
-0.878931188971585, 0.665670601363721, -0.665670601363721, -0.787856768435423,
0.787856768435423), 2, 3),
residuals = matrix(c(0.0367718470212829,
-0.0367718470212828, -0.214589957886459, 0.214589957886458, -0.647955514235535,
0.647955514235536), 2, 3,
dimnames = list(t=2:3, i=1:3)),
df.residual = 2*3 - 3,
logLik = 0.704787094195681,
design = "balanced",
call = quote(opm(x = y~x1+x2, data = d, n.samp = 10)),
.Environment = environment(),
time.indicators = FALSE,
index = c('i', 't'),
terms = structure(y~x1+x2,
variables = quote(list(y, x1, x2)),
factors = matrix(c(0, 1, 0, 0, 0, 1), 3, 2,
dimnames=list(c('y', 'x1', 'x2'),
c('x1', 'x2'))),
term.labels = c('x1', 'x2'),
order = c(1L, 1L),
intercept = 0L,
response = 1L,
.Environment = parent.frame(),
predvars = quote(list(y, x1, x2)),
dataClasses = c(y='numeric', x1='numeric', x2='numeric'),
class=c('terms', 'formula'))),
class = 'opm')
## default index is the first two columns of the data
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
expect_equal(opm(y~x1+x2, d,
n.samp = 10),
expected)
d <- d[,c('y', 'x1', 't', 'i', 'x2')]
## specify index column by position
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
expected$call <- quote(opm(x = y~x1+x2, data = d,
index = 4:3, n.samp = 10))
expect_equal(opm(y~x1+x2, d, index=4:3,
n.samp = 10),
expected)
## specify index column by name
set.seed(123)
expected$call <- quote(opm(x = y~x1+x2, data = d,
index = c('i', 't'), n.samp = 10))
expect_equal(opm(y~x1+x2, d, index=c('i', 't'), n.samp = 10),
expected)
})
test_that('early dropouts', {
suppressWarnings(RNGversion("3.5.0"))
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion")
xx <- array(NA, dim=dim(x) + c(0, 0, 1))
xx[,,seq_len(ncol(y))] <- x
yy <- matrix(NA, nrow(y), ncol(y)+1)
yy[,seq_len(ncol(y))] <- y
## one more case, but immediately drops out
expect_equal(dim(xx), dim(x)+c(0,0,1))
expect_equal(dim(yy), dim(y)+c(0,1))
expect_equal(opm(xx, yy, n = 10),
structure(list(samples = list(rho = c(0.805, 0.425,
-0.181, 0.763,
0.581, 0.965,
0.0559999999999999, 0.533,
0.648, -0.086),
sig2 = 1 / c(11.6175662271191,
0.919115793520271,
0.989369792715488,
0.00440918521157215,
1.35621264272223,
4.24440799476086,
0.0931249250470442,
1.56847551750474,
0.371922970545138,
5.00022577787265),
beta = matrix(c(0.149966256134789, -1.07201170366086,
-3.32784640301039, 6.31976053552415,
0.399370400131068, -0.409833517459282,
0.942409502841064, -1.01613441888489,
-1.06148590924992, -1.70364397749032,
-0.746392848950218, -1.22143920569537,
-1.63866906134468, 7.01947247537598,
-0.536929048279301, -0.845295844019301,
-0.565341215754439, 0.0765501718682319,
-0.171463431207414, -1.13494915857112),
10, 2)),
fitted.values = matrix(c(0.878931188971585, -0.878931188971585,
0.665670601363721, -0.665670601363721,
-0.787856768435423, 0.787856768435423,
0, 0), 2, 4),
residuals = matrix(c(0.0367718470212829, -0.0367718470212828,
-0.214589957886459, 0.214589957886458,
-0.647955514235535, 0.647955514235536,
0, 0), 2, 4),
df.residual = 2*3 - 3,
logLik = 0.637475935193608,
design = "unbalanced (with dropouts)",
call = quote(opm(x = xx, y = yy, n.samp = 10)),
.Environment = environment(),
time.indicators = FALSE),
class = 'opm'))
})
test_that('confint', {
expect_equal(confint(o),
matrix(c(-0.159625, 0.111707292991574,
-2.96240085726837, -1.54479234382358,
0.929, 178.18554581309,
5.10985655317045, 5.45731495708674),
nrow = 4, ncol = 2,
dimnames = list(c("rho", "sig2", "beta.x1", "beta.x2"),
c("2.5%", "97.5%"))))
expect_equal(confint(o, 1:2, level = 0.68),
matrix(c(-0.0235200000000001, 0.215660754512765,
0.78652, 7.19646832807839),
nrow = 2, ncol = 2,
dimnames = list(c("rho", "sig2"), c("16%", "84%"))))
})
test_that('coef', {
expect_equal(coef(o),
c(rho = 0.557,
sig2 = 0.874045960780248,
beta.x1 = -0.712983968172086,
beta.x2 = -0.655867032352328))
expect_error(coef(o, probs = .6), "Arguments 'probs' and 'names' are not allowed")
expect_error(coef(o, names = TRUE), "Arguments 'probs' and 'names' are not allowed")
})
test_that('print', {
expect_equal(capture.output(o),
c('\tPanel design: balanced',
'',
'Call:',
'opm(x = y ~ x1 + x2, data = d, n.samp = 10)',
'',
'Coefficients:',
' mean (SD) med 95-CI',
'rho 0.450900 (0.39) 0.55700 (-0.16, 0.93)',
'sig2 24.422159 (71.18) 0.87405 (0.11, 178.19)',
'beta.x1 -0.077945 (2.55) -0.71298 (-3.0, 5.1)',
'beta.x2 0.023554 (2.51) -0.65587 (-1.5, 5.5)',
''))
})
test_that('summary', {
expect_equal(capture.output(summary(o)),
c('Call:',
'opm(x = y ~ x1 + x2, data = d, n.samp = 10)',
'',
'Parameter estimates:',
' <--95CI <--68CI med 68CI--> 95CI-->',
'rho -0.15962 -0.02352 0.55700 0.786520 0.9290',
'sig2 0.11171 0.21566 0.87405 7.196468 178.1855',
'beta.x1 -2.96240 -1.42573 -0.71298 0.703472 5.1099',
'beta.x2 -1.54479 -1.18338 -0.65587 -0.032576 5.4573'))
})
test_that('DIC', {
expect_equal(DIC(o),
24.8318940856621)
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.