Nothing
tests/testthat/test_models.r
In msm: Multi-State Markov and Hidden Markov Models in Continuous Time
context("msm simple model likelihoods")
test_that("simple model, exact death times",{
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
method="BFGS", control=list(trace=5, REPORT=1))
print(cav.msm)
printold.msm(cav.msm)
expect_equal(4908.81676837903, cav.msm$minus2loglik)
})
test_that("simple model, not exact death times",{
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = FALSE, fixedpars=TRUE)
expect_equal(4833.00640639644, cav.msm$minus2loglik)
})
test_that("simple model, covariates",{
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
covariates = ~ sex, covinits = list(sex=rep(0.01, 7)))
expect_equal(4909.17147259115, cav.msm$minus2loglik)
})
test_that("autogenerated inits reproduce crudeinits",{
cinits <- crudeinits.msm(state ~ years, PTNUM, data=cav, qmatrix=twoway4.q)
expect_equal(msm( state ~ years, subject=PTNUM, data = cav, qmatrix = cinits, deathexact = TRUE, fixedpars=TRUE)$minus2loglik,
msm( state ~ years, subject=PTNUM, data = cav, qmatrix=twoway4.i, gen.inits=TRUE, deathexact = TRUE, fixedpars=TRUE)$minus2loglik)
})
test_that("data as global variables",{
state.g <- cav$state; time.g <- cav$years; subj.g <- cav$PTNUM
cav.msm <- msm(state.g ~ time.g, subject=subj.g, qmatrix = twoway4.i, gen.inits=TRUE, fixedpars=TRUE)
expect_equal(4119.9736299032, cav.msm$minus2loglik)
})
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), control=list(fnscale=1))
test_that("psor model: covariates, constraints",{
print(psor.msm)
printold.msm(psor.msm)
expect_equal(1114.89946121717, psor.msm$minus2loglik, tol=1e-06)
expect_equal(0.0959350004999946, psor.msm$Qmatrices$baseline[1,2], tol=1e-06)
expect_equal(exp(psor.msm$Qmatrices$logbaseline[c(5,10,15)]), psor.msm$Qmatrices$baseline[c(5,10,15)], tol=1e-06)
})
test_that("psor model: transition-specific covariates",{
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, covariates = list("2-3"=~ollwsdrt, "3-4"=~hieffusn), control=list(fnscale=1))
expect_equal(hazard.msm(psor.msm)$ollwsdrt["State 1 - State 2","HR"], 1)
expect_equal(hazard.msm(psor.msm)$hieffusn["State 2 - State 3","HR"], 1)
expect_error(msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = list("meep"=~ollwsdrt, "3-4"=~hieffusn), control=list(fnscale=1)),
"not in format \"number-number\"")
})
test_that("psor model: transition-specific covariates with pci",{
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = list("2-3"=~ollwsdrt, "3-4"=~hieffusn),
pci = 10, fixedpars=7:8, control=list(fnscale=1))
h <- hazard.msm(psor.msm)
expect_equal(h$ollwsdrt["State 1 - State 2","HR"], 1)
expect_equal(h$hieffusn["State 1 - State 2","HR"], 1)
expect_true(!isTRUE(all.equal(h$"timeperiod[10,Inf)"["State 1 - State 2","HR"], 1)))
expect_equal(h$"timeperiod[10,Inf)"["State 2 - State 3","HR"], 1)
expect_equal(h$"timeperiod[10,Inf)"["State 3 - State 4","HR"], 1)
})
psor.nocen.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), center=FALSE)
test_that("no covariate centering",{
expect_equal(exp(psor.nocen.msm$Qmatrices$logbaseline[c(5,10,15)]), psor.nocen.msm$Qmatrices$baseline[c(5,10,15)], tol=1e-06)
})
test_that("gen.inits with missing values for state / time",{
psor2 <- psor; psor2$ptnum[13:14] <- psor2$months[7:8] <- psor2$state[7:8] <- NA
psor2.msm <- msm(state ~ months, subject=ptnum, data=psor2, gen.inits=TRUE, qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), fixedpars=TRUE)
expect_equal(1179.95441284169, psor2.msm$minus2loglik, tol=1e-06)
})
test_that("exact transition times using exacttimes",{
msmtest5 <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, exacttimes=TRUE, fixedpars=TRUE)
expect_equal(3057.85781916437, msmtest5$minus2loglik, tol=1e-06)
})
test_that("exact transition times using obstype vector of all 2",{
msmtest5 <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, obstype=rep(2, nrow(bos)), fixedpars=TRUE)
expect_equal(3057.85781916437, msmtest5$minus2loglik, tol=1e-06)
})
test_that("exact transition times with death, should be same",{
expect_warning(msmtested <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, deathexact=5, obstype=rep(2, nrow(bos)), exacttimes=TRUE, fixedpars=TRUE), "Ignoring death")
expect_equal(3057.85781916437, msmtested$minus2loglik, tol=1e-06)
(msmtest5 <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, deathexact=5, obstype=rep(2, nrow(bos)), fixedpars=TRUE)) # no warning, inconsistently
expect_equal(3057.85781916437, msmtest5$minus2loglik, tol=1e-06)
})
cav$statefac <- factor(cav$state)
cav$statefac2 <- factor(cav$state, labels=c("none","mild","severe","death"))
test_that("factors as states, death",{
cav.msm <- msm( statefac ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
method="BFGS", control=list(trace=5, REPORT=1))
expect_equal(4908.81676837903, cav.msm$minus2loglik)
expect_error(msm( statefac2 ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "state variable should be numeric or a factor with ordinal numbers as levels")
})
test_that("factor covariates with factor() in formula",{
## Should be no need for users to do this. factors should already be identified as such in the data frame
expect_warning(cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ factor(pdiag), covinits=list(sex=rep(0.1,7)), fixedpars=TRUE), "covariate .+ unknown")
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ factor(pdiag), covinits=list("factor(pdiag)Hyper"=rep(0.1,7)), fixedpars=TRUE)
expect_equal(4793.20440566637, cavfaccov.msm$minus2loglik, tol=1e-06)
})
test_that("factors as global variables",{
state.g <- cav$state; time.g <- cav$years; subj.g <- cav$PTNUM; pdiag.g <- factor(cav$pdiag)
cavfaccov.msm <- msm(state.g ~ time.g, subject=subj.g, qmatrix = twoway4.q, covariates = ~ pdiag.g, fixedpars=TRUE)
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
})
test_that("factor covariates using existing factors: inits are given to contrasts ",{
## Warnings could be more informative here
expect_warning(cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ pdiag, covinits=list(pdiag=rep(0.1,7)), fixedpars=TRUE), "covariate .+ unknown")
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
expect_warning(cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ pdiag, covinits=list(pdiagNonexistentlevel=rep(0.1,7)), fixedpars=TRUE), "covariate .+ unknown")
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ pdiag, covinits=list(pdiagHyper=rep(0.1,7)), fixedpars=TRUE)
expect_equal(4793.20440566637, cavfaccov.msm$minus2loglik, tol=1e-06)
cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ pdiag,
covinits=list(pdiagHyper=rep(0.1,7),pdiagIDC=rep(0.1,7),pdiagIHD=rep(0.1,7),pdiagOther=rep(0.1,7),pdiagRestr=rep(0.1,7)), fixedpars=TRUE) # OK
expect_equal(4793.13035886505, cavfaccov.msm$minus2loglik, tol=1e-06)
})
context("censored states")
test_that("censored states: final state censored",{
cavcens.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens, qmatrix=twoway4.q, censor=99, fixedpars=TRUE)
expect_equal(4724.26606344485, cavcens.msm$minus2loglik, tol=1e-06)
v <- viterbi.msm(cavcens.msm)
expect_equal(v$observed[v$observed<10], v$fitted[v$observed<10])
expect_equal(v$fitted[v$observed==99][1], 3)
})
test_that("two kinds of censoring",{
cavcens2.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens2, qmatrix=twoway4.q, censor=c(99, 999), censor.states=list(c(1,2,3), c(2,3)), fixedpars=TRUE)
expect_equal(4678.23348518727, cavcens2.msm$minus2loglik, tol=1e-06)
v <- viterbi.msm(cavcens2.msm)
expect_equal(v$observed[v$observed<10], v$fitted[v$observed<10])
expect_equal(v$fitted[v$observed==99][1], 3)
})
test_that("intermediate state censored",{
cavcens3.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens3, qmatrix=twoway4.q, censor=c(99, 999), censor.states=list(c(2,3), c(1,2,3)), fixedpars=TRUE)
expect_equal(4680.66073438518, cavcens3.msm$minus2loglik, tol=1e-06)
v <- viterbi.msm(cavcens3.msm)
expect_equal(v$observed[v$observed<10], v$fitted[v$observed<10])
expect_true(all(v$fitted[v$observed==99] %in% 2:3))
expect_true(all(v$fitted[v$observed==999] %in% 1:3))
})
test_that("first state censored",{
cav.cens4 <- cav
cav.cens4$state[c(1,8,12,22)] <- 99
cavcens4.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens4, qmatrix=twoway4.q,
censor=c(99), censor.states=list(c(2,3)), fixedpars=TRUE)
expect_equal(4846.06045097812, cavcens4.msm$minus2loglik)
v <- viterbi.msm(cavcens4.msm)
expect_true(all(v$fitted[v$observed==99] %in% 2:3))
})
test_that("piecewise constant intensities with pci",{
cav5.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
pci = c(5), covinits = list("timeperiod[5,Inf)"=rep(0.01,7)),
)
expect_equal(4906.01423796805, cav5.msm$minus2loglik, tol=1e-06)
cav10.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, pci = c(5,10), fixedpars=TRUE,
covinits = list("timeperiod[5,10)"=rep(0.01,7), "timeperiod[10,Inf)"=rep(0.01,7)),
)
expect_equal(4905.61646158639, cav10.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with pci, cut points outside data",{
## Make sure works for pci outside time range, with warning
expect_warning(cav5.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE, pci = c(-1,5,50,60), covinits = list("timeperiod[5,Inf)"=rep(0.01,7))), "cut point .+ less than")
expect_warning(cav.pci.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE, pci = c(-1,50,60)))
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE)
expect_equal(cav.pci.msm$minus2loglik, cav.msm$minus2loglik) # degrades to time-homogeneous if all cuts outside data
})
suppressWarnings(RNGversion("3.5.0"))
set.seed(22061976)
cav$pdiag3 <- cav$pdiag
cav$pdiag3[!cav$pdiag %in% c("IDC","IHD")] <- "Other"
cav$pdiag3 <- factor(cav$pdiag3)
subs <- cav$PTNUM %in% sample(unique(cav$PTNUM), 50)
cavsub <- subset(cav, subs)
cavsub$maxtime <- tapply(cavsub$years, cavsub$PTNUM, max)[as.character(cavsub$PTNUM)]
cavsub.extra <- cavsub[cavsub$years==0 & cavsub$maxtime >= 5,]
cavsub.extra$years <- 5
cavsub.extra$state <- 99
cavsub.extra10 <- cavsub[cavsub$years==0 & cavsub$maxtime >= 10,]
cavsub.extra10$years <- 10
cavsub.extra10$state <- 999
cavsub2 <- rbind(cavsub, cavsub.extra, cavsub.extra10)
cavsub2 <- cavsub2[order(cavsub2$PTNUM, cavsub2$years),]
cavsub2$after5 <- ifelse(cavsub2$years>=5 & cavsub2$years<10, 1, 0)
cavsub2$after10 <- ifelse(cavsub2$years>=10, 1, 0)
cavsub2$after510 <- as.numeric(cavsub2$after5 | cavsub2$after10)
test_that("piecewise constant intensities with pci, with other covariates",{
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub,
qmatrix = twoway4.q, covariates = ~ pdiag3 + sex, deathexact = TRUE, pci = c(5,10), fixedpars=TRUE,
covinits = list("timeperiod[5,10)"=rep(0.01,7), "timeperiod[10,Inf)"=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
)
expect_equal(448.122802051545, cav5cov.msm$minus2loglik, tol=1e-06)
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub2,
qmatrix = twoway4.q, deathexact = TRUE, covariates = ~ pdiag3 + sex + after5 + after10,
censor=c(99,999), censor.states=list(1:4,1:4),
covinits = list(after5=rep(0.01,7), after10=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
fixedpars=TRUE
)
expect_equal(448.122802051545, cav5cov.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with pci, with uncentered covariates",{
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub,
qmatrix = twoway4.q, covariates = ~ pdiag3 + sex, deathexact = TRUE, pci = c(5,10), fixedpars=TRUE, center=FALSE,
covinits = list("timeperiod[5,10)"=rep(0.01,7), "timeperiod[10,Inf)"=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
)
expect_equal(449.691983558378, cav5cov.msm$minus2loglik, tol=1e-06)
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub2,
qmatrix = twoway4.q, deathexact = TRUE, covariates = ~ pdiag3 + sex + after5 + after10,
censor=c(99,999), censor.states=list(1:4,1:4), center=FALSE,
covinits = list(after5=rep(0.01,7), after10=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
fixedpars=TRUE
)
expect_equal(449.691983558378, cav5cov.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with pci, with other censored states",{
cav.cens <- cav
cav.cens$state[cav.cens$state==4][1:50] <- 99
cav5cens.msm <- msm(state ~ years, subject=PTNUM, data = cav.cens,
qmatrix = twoway4.q, deathexact = TRUE, censor=99, pci = 5,
covinits=list("timeperiod[5,Inf)"=rep(0.02,7)), fixedpars=TRUE,
method="BFGS", control=list(trace=5, REPORT=1))
expect_equal(4754.41981265159, cav5cens.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with covariates in HMMs",{
misccov.msm <- msm(state ~ years, subject = PTNUM, data = cav,
qmatrix = oneway4.q, ematrix=ematrix, deathexact = 4, fixedpars=TRUE,
pci = 5, covinits=list("timeperiod[5,Inf)"=rep(0.0001,5)),
misccovariates = ~dage + sex, misccovinits = list(dage=c(0.01,0.02,0.03,0.04), sex=c(-0.013,-0.014,-0.015,-0.016)))
expect_equal(4306.29865288466, misccov.msm$minus2loglik, tol=1e-06)
})
context("output functions")
test_that("qmatrix.msm for psor model, defaults",{
expect_equal(c(-0.0959350004999946, 0, 0, 0, 0.0959350004999946, -0.164306508892574, 0, 0, 0, 0.164306508892574, -0.254382807485639, 0, 0, 0, 0.254382807485639, 0), as.numeric(qmatrix.msm(psor.msm)$estimates), tol=1e-03)
expect_equal(c(0.0115942726096754, 0, 0, 0, 0.0115942726096754, 0.0196169975000406, 0, 0, 0, 0.0196169975000406, 0.0375066077515386, 0, 0, 0, 0.0375066077515386, 0), as.numeric(qmatrix.msm(psor.msm)$SE), tol=1e-03)
expect_error(qmatrix.msm(psor.msm, ci="normal", cl=-1), "expected cl")
expect_equivalent(qmatrix.msm(psor.msm, sojourn=TRUE)$sojourn[1:3],
sojourn.msm(psor.msm)$estimates)
qmatrix.msm(psor.msm, ci="normal", B=2)
qmatrix.msm(psor.msm, sojourn=TRUE, ci="normal", B=2)
expect_error(qmatrix.msm("foo"), "expected .+ msm model")
expect_error(qmatrix.msm(psor.msm, covariates="foo"), "covariates argument must be")
expect_warning(qmatrix.msm(psor.msm, covariates=list(foo=1)), "ignoring")
})
test_that("qmatrix.msm defaults to Qmatrices in object",{
expect_equal(psor.msm$Qmatrices$baseline, qmatrix.msm(psor.msm, covariates="mean", ci="none"))
expect_equal(psor.nocen.msm$Qmatrices$baseline, qmatrix.msm(psor.nocen.msm, covariates=0, ci="none"))
})
test_that("qmatrix.msm with supplied covariates",{
qmat <- qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4))
expect_equal(c(-0.121430585652200, 0, 0, 0, 0.121430585652200, -0.207972362475868, 0, 0, 0, 0.207972362475868, -0.257535341208494, 0, 0, 0, 0.257535341208494, 0), as.numeric(qmat$estimates), tol=1e-03)
expect_equal(c(0.0162156605802465, 0, 0, 0, 0.0162156605802465, 0.0266727053124233, 0, 0, 0, 0.0266727053124233, 0.0364321127089265, 0, 0, 0, 0.0364321127089265, 0), as.numeric(qmat$SE), tol=1e-04)
expect_warning(qmatrix.msm(psor.msm, covariates=list(hieffusn=0.1, foo=0.4)), "Covariate .+ unknown")
})
test_that("qmatrix.msm with non-default confidence limits",{
qmat <- qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4), cl=0.99)
expect_equal(c(-0.171282667596986, 0, 0, 0, 0.0860880282792585, -0.289385121267802, 0, 0, 0, 0.149463467106753, -0.370756460718086, 0, 0, 0, 0.178889538008097, 0), as.numeric(qmat$L), tol=1e-04)
})
test_that("qmatrix.msm sojourn component",{
soj <- qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4), sojourn=TRUE)$sojourn
expect_equal(c(8.23515751512713, 4.80833120370037, 3.88296221911705, Inf), as.numeric(soj), tol=1e-03)
expect_equal(as.numeric(soj[1:3]), sojourn.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4))[,"estimates"])
})
test_that("qmatrix.msm bug for user-supplied covariates fixed in 1.1.3",{
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$SE[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0, ollwsdrt=0))$SE[1,2])
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$SE[1,1],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0, ollwsdrt=0))$SE[1,1])
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$L[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0, ollwsdrt=0))$L[1,2])
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$L[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0))$L[1,2])
cm <- psor.nocen.msm$qcmodel$covmeans
expect_equal(qmatrix.msm(psor.nocen.msm, covariates="mean")$SE[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=cm["hieffusn"], ollwsdrt=cm["ollwsdrt"]))$SE[1,2])
})
test_that("qmatrix.msm: unspecified covariate values default to zero",{
expect_equal(psor.nocen.msm$Qmatrices$baseline, qmatrix.msm(psor.nocen.msm, covariates=list(ollwsdrt=0), ci="none")) # missing covs default to zero
expect_equal(qmatrix.msm(psor.msm, covariates=list(hieffusn=0)),
qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0, hieffusn=0)))
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0)),
qmatrix.msm(psor.nocen.msm, covariates=list(ollwsdrt=0, hieffusn=0)))
})
test_that("sojourn.msm",{
expect_equal(psor.msm$sojourn, sojourn.msm(psor.msm, covariates="mean"))
expect_equal(psor.nocen.msm$sojourn, sojourn.msm(psor.nocen.msm, covariates=0))
soj <- sojourn.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4))
expect_equal(c(8.23515751512713, 4.80833120370037, 3.88296221911705, 1.09971073904434, 0.616674252838334, 0.549301375677405, 6.33875136203292, 3.73961380505919, 2.94271599303942, 10.6989240349703, 6.18246967994404, 5.12363260020806), as.numeric(unlist(soj)), tol=1e-04)
soj <- sojourn.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4), cl=0.99)
expect_equal(5.83830234564607, soj[1,"L"], tol=1e-04)
expect_error(sojourn.msm("foo"), "expected .+ msm model")
expect_error(sojourn.msm(psor.msm, covariates="foo"),"covariates argument must be")
expect_warning(sojourn.msm(psor.msm, covariates=list(foo=1)), "ignoring")
})
test_that("pmatrix.msm",{
expect_equal(0.149287738928777, pmatrix.msm(psor.msm, ci="none", t=10)[1,3], tol=1e-04)
p <- pmatrix.msm(psor.msm, t=10, covariates=list(ollwsdrt=0.1, hieffusn=0.2))
expect_equal(0.18196160265907, p[1,3], tol=1e-04)
set.seed(22061976); expect_equal(0.12, pmatrix.msm(psor.msm, ci="normal", B=3)$L[2,3], tol=1e-01, scale=1)
expect_error(pmatrix.msm("foo"), "expected .+ msm model")
expect_error(pmatrix.msm(psor.msm, t="foo"), "must be a positive number")
expect_error(pmatrix.msm(psor.msm, -9), "must be a positive number")
expect_equivalent(unclass(pmatrix.msm(psor.msm, 0)), diag(4))
expect_warning(pmatrix.msm(psor.msm, 1, covariates=list(foo=1)), "ignoring")
})
test_that("pnext.msm",{
expect_equal(pnext.msm(psor.msm, ci="none")$estimate[1,2], 1)
expect_equal(pnext.msm(psor.msm, ci="delta")$estimate[1,2], 1)
expect_equal(pnext.msm(psor.msm, ci="normal", B=2)$estimate[1,2], 1)
})
test_that("qratio.msm",{
q <- qratio.msm(psor.msm, c(1,2), c(2,3))
expect_equal(c(0.583878474075081, 0.0996029045389022, 0.417943274168735, 0.815694601537263), as.numeric(q), tol=1e-04)
q <- qratio.msm(psor.msm, c(1,2), c(2,3), cl=0.99)
expect_equal(0.376262194364283, as.numeric(q["L"]), tol=1e-04)
q <- qratio.msm(psor.msm, c(1,1), c(2,3))
expect_equal(c(-0.583878474075081, 0.0996029045389022, -0.815694601537263, -0.417943274168735), as.numeric(q), tol=1e-04)
q <- qratio.msm(psor.msm, c(2,2), c(2,3))
expect_equivalent(c(-1,0,-1,-1), q)
qratio.msm(psor.msm, c(1,2), c(2,3), ci="norm", B=2)
expect_error(qratio.msm("foo"))
expect_error(qratio.msm(psor.msm, "foo"))
expect_error(qratio.msm(psor.msm, c(1,8), c(1,0)))
expect_error(qratio.msm(psor.msm, c(1,2), c(1,0)))
expect_error(qratio.msm(psor.msm, c(1,2), c(2,3), cl="foo"))
expect_error(qratio.msm(psor.msm, c(1,2), c(2,3), cl=2), "expected cl in")
qq <- qratio.msm(psor.msm, c(1,2), c(2,2))
QQ <- qmatrix.msm(psor.msm)
expect_equal(qq[["estimate"]], (QQ[1,2]/QQ[2,2])[["estimate"]])
qq <- qratio.msm(psor.msm, c(1,1), c(2,2))
QQ <- qmatrix.msm(psor.msm)
expect_equal(qq[["estimate"]], (QQ[1,1]/QQ[2,2])[["estimate"]])
})
test_that("coef.msm",{
co <- coef.msm(psor.msm)
expect_equal(0.498319866154661, co$hieffusn[1,2], tol=1e-04)
expect_error(coef.msm("foo"))
})
test_that("hazard.msm",{
haz <- hazard.msm(psor.msm)
expect_equal(0.385347226135311, haz$ollwsdrt[1,2], tol=1e-04)
expect_equal(0.385347226135311, haz$ollwsdrt[2,2], tol=1e-04)
expect_equal(2.35928404626333, haz$hieffusn[1,3], tol=1e-04)
expect_equal(2.35928404626333, haz$hieffusn[3,3], tol=1e-04)
haz <- hazard.msm(psor.msm, hazard.scale=2)
expect_equal(0.148492484690178, haz$ollwsdrt[1,2], tol=1e-04)
expect_equal(haz$ollwsdrt[1,2], haz$ollwsdrt[2,2])
expect_equal(haz$hieffusn[1,3], haz$hieffusn[3,3])
haz <- hazard.msm(psor.msm, hazard.scale=c(1,2))
expect_equal(0.385347226135311, haz$ollwsdrt[1,2], tol=1e-04)
expect_equal(haz$ollwsdrt[1,2], haz$ollwsdrt[2,2], tol=1e-04)
expect_equal(haz$hieffusn[1,3], haz$hieffusn[3,3])
expect_error(hazard.msm("foo"))
expect_error(hazard.msm(psor.msm, hazard.scale="foo"))
expect_error(hazard.msm(psor.msm, hazard.scale=c(1,2,3)), "hazard.scale of length")
})
test_that("transient.msm and absorbing.msm",{
expect_equivalent(c(1,2,3), transient.msm(psor.msm))
expect_equivalent(4, absorbing.msm(psor.msm))
expect_error(transient.msm("foo"))
expect_error(absorbing.msm("foo"))
expect_error(transient.msm(qmatrix="foo"))
expect_error(absorbing.msm(qmatrix="foo"))
expect_error(transient.msm(qmatrix=c(1,4,5,6)))
expect_error(absorbing.msm(qmatrix=c(1,4,5,6)))
expect_error(transient.msm(qmatrix=cbind(c(1,2,3),c(1,3,2))))
expect_error(absorbing.msm(qmatrix=cbind(c(1,2,3),c(1,3,2))))
expect_error(transient.msm())
expect_error(absorbing.msm())
})
test_that("prevalence.msm",{
p <- prevalence.msm(psor.msm)
expect_equal(59, p$Observed[5,5], tol=1e-06)
expect_equal(59, p$Expected[5,5], tol=1e-06)
expect_equal(57.35294, p$"Observed percentages"[4,4], tol=1e-03)
expect_equal(49.96882, p$"Expected percentages"[4,4], tol=1e-03)
summ <- summary.msm(psor.msm)
expect_equal(summ$prevalences, p)
p <- prevalence.msm(psor.msm, times=seq(0,60,5))
expect_equal(63, p$Observed[5,5], tol=1e-06)
expect_equal(63, p$Expected[5,5], tol=1e-06)
expect_equal(50.70423, p$"Observed percentages"[4,4], tol=1e-03)
expect_equal(41.46338, p$"Expected percentages"[4,4], tol=1e-03)
expect_error(prevalence.msm("foo"))
expect_error(summary.msm("foo"))
p <- prevalence.msm(psor.msm, covariates=list(hieffusn=0, ollwsdrt=1))
expect_equal(p$Observed[1,1], 1)
p <- prevalence.msm(psor.msm, covariates=list(hieffusn=0, ollwsdrt=1),
ci="normal", B=10)
expect_true(is.numeric(p[[2]]$ci[,,"2.5%"][1,1]))
plot.prevalence.msm(psor.msm)
})
test_that("pmatrix.piecewise.msm",{
times <- c(5, 10, 15)
covariates <- list(list(ollwsdrt=0, hieffusn=0),
list(ollwsdrt=0, hieffusn=1),
list(ollwsdrt=1, hieffusn=0),
list(ollwsdrt=1, hieffusn=1)
)
p <- pmatrix.msm(psor.msm, 3, covariates=covariates[[1]])
pp <- pmatrix.piecewise.msm(psor.msm, 0, 3, times, covariates)
expect_equal(pp[1,3], p[1,3], tol=1e-04)
p <- pmatrix.piecewise.msm(psor.msm, 0, 7, times, covariates)
expect_equal(0.172773087945103, p[1,3], tol=1e-04)
pp <- pmatrix.piecewise.msm(psor.msm, 0, 19, times, covariates)
expect_equal(0.0510873669808412, pp[1,3], tol=1e-04)
p <- pmatrix.msm(psor.msm, 5, covariates=covariates[[1]]) %*% pmatrix.msm(psor.msm, 5, covariates=covariates[[2]]) %*% pmatrix.msm(psor.msm, 5, covariates=covariates[[3]]) %*% pmatrix.msm(psor.msm, 4, covariates=covariates[[4]])
expect_equal(pp[1,3], p[1,3], tol=1e-04)
covariates <- list(list(ollwsdrt=0, hieffusn=0),list(ollwsdrt=0, hieffusn=1))
p <- pmatrix.piecewise.msm(psor.msm, 0, 7, times=5, covariates)
expect_equal(0.172773087945103, p[1,3], tol=1e-04)
expect_error(pmatrix.piecewise.msm("foo", 1,2, c(1, 2), c(1,2)), "expected .+ msm model")
expect_error(pmatrix.piecewise.msm("foo", 1, 2, c(1, 0.5, 2), list(0, 1, 0, 1)), "expected .+ msm model")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 0.5, 2), list(0, 1, 0, 1)), "times should be a vector of numbers in increasing order")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, "rubbish", list(0, 1, 0, 1)),"times should be a vector of numbers in increasing order")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 1.5, 2), "rubbish"))
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 1.5, 2), list("rubbish", "foo","bar","boing")), "covariates argument")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 1.5, 2), list(0, 1, 0, 1)), "covariates argument")
})
test_that("pmatrix.piecewise.msm given just Q matrices",{
p1 <- pmatrix.piecewise.msm(
t1 = 0,
t2 = 6,
times = c(2, 4),
covariates = list(list(cov1 = 1), list(cov1 = 2), list(cov1 = 3)),
qlist = list(
Q1 = rbind(c(0.9, 0.1), c(0.1, 0.9)),
Q2 = rbind(c(0.8, 0.2), c(0.1, 0.9)),
Q3 = rbind(c(0.7, 0.3), c(0.1, 0.9))
)
)
expect_equal(round(p1), rbind(c(166, 238), c(99, 304)))
})
test_that("totlos.msm",{
tl <- totlos.msm(psor.msm, ci="normal", B=2)
sj <- sojourn.msm(psor.msm, ci="none")
expect_equal(tl[1,1], sj$estimates[[1]])
en <- envisits.msm(psor.msm, ci="normal", B=2)
expect_equal(en[1,2], 1)
})
test_that("logLik.msm",{
expect_equivalent(unclass(logLik.msm(psor.msm)), psor.msm$minus2loglik / -2)
expect_error(logLik.msm("foo"))
})
test_that("lrtest.msm",{
psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~ollwsdrt, constraint = list(ollwsdrt=c(1,1,2)), control=list(fnscale=1))
expect_equal(lrtest.msm(psor2.msm, psor.msm)[1,"-2 log LR"],
psor2.msm$minus2loglik - psor.msm$minus2loglik)
})
test_that("qmatrix subset function",{
Q <- qmatrix.msm(psor.msm)
expect_equivalent(Q[1,2]["SE"], Q$SE[1,2])
expect_error(Q[1,2,3,4,5], "unused arguments")
expect_error(Q[1], "Two dimensions must be supplied")
expect_equivalent(Q[], Q)
Q[1,]
Q[,2]
Q[c(1,2),]
Q[c(2,1),]
Q[c(1,2),c(1,2)]
})
test_that("efpt.msm",{
Q <- twoway4.q
expect_equal(efpt.msm(qmatrix=Q, tostate=3), c(Inf,Inf,0,Inf))
Q <- rbind(c(-0.25,0.25,0), c(0.166, -0.332, 0.166), c(0, 0.25, -0.25))
expect_equal(efpt.msm(qmatrix=Q, tostate=3)[c(1,2)], solve(-Q[1:2,1:2], c(1,1)),tol=1e-06)
Q <- twoway4.q; Q[2,4] <- Q[2,1] <- 0; diag(Q) <- 0; diag(Q) <- -rowSums(Q)
expect_equal(efpt.msm(qmatrix=Q, tostate=3), c(Inf, 6.02409638554217, 0, Inf), tol=1e-05)
expect_equal(efpt.msm(psor.msm, tostate=c(2)), c(10.4237243422138, 0, Inf, Inf), tol=1e-05)
expect_equal(efpt.msm(psor.msm, tostate=c(3)), c(16.5099104995137, 6.08618615729989, 0, Inf), tol=1e-05)
expect_equal(efpt.msm(psor.msm, tostate=c(2,3)), c(10.4237243422138, 0, 0, Inf), tol=1e-05)
})
test_that("ppass.msm",{
pp <- ppass.msm(psor.msm, tot=10)
pm <- pmatrix.msm(psor.msm, t=10)
expect_equal(pp[,4], pm[,4]) # state 4 is absorbing
pp <- ppass.msm(qmatrix=twoway4.q, tot=1000)
expect_equal(pp[1,2], 0.5)
expect_warning(ppass.msm(qmatrix=twoway4.q, tot=100, ci="normal"), "No fitted model supplied: not calculating confidence intervals")
})
test_that("score residuals",{
sres <- scoreresid.msm(psor.msm)
expect_equal(c(0.0608163009112361, 0.187998750251689, 0.0143302186951471), as.numeric(sres[1:3]), tol=1e-05)
psor2 <- na.omit(psor); psor2$months[psor2$ptnum==5] <- psor2$months[psor2$ptnum==5]*10
psor.0.q <- rbind(c(0,0.1,0,0),c(0,0,0.2,0),c(0,0,0,0.3),c(0,0,0,0))
psor.infl.msm <- msm(state ~ months, subject=ptnum, data=psor2, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), qmatrix = psor.0.q)
sres <- scoreresid.msm(psor.infl.msm)
if (interactive()) sres <- scoreresid.msm(psor.infl.msm, plot=TRUE)
expect_equal(names(which.max(sres)), "5")
})
test_that("observed, expected etc",{
## two covariates
expect_equal(get.covhist(psor.msm)$example$time, c(6.4606, 17.078, 26.3217, 30.5763, 6.4052, 7.6893, 3.3593, 12.7775))
expect_equivalent(observed.msm(psor.msm)$obstab[,"State 4"], c(0, 6, 25, 39, 48, 50, 52, 53, 55, 56, 57))
expect_equal(as.numeric(expected.msm(psor.msm, covariates="population")$Expected[1:5,"State 3"]),
c(0, 11.7058415609289, 14.8078062584701, 9.50355791490463, 5.83372436779321), tol=1e-05)
expect_equal(as.numeric(expected.msm(psor.msm, covariates="mean")$Expected[1:5,"State 3"]),
c(0, 11.165395364321, 14.8732589026358, 9.67192721319246, 6.26754773418822), tol=1e-05)
})
test_that("observed, expected etc with one covariate, should run",{
expect_error({
psor1.msm <- msm(state ~ months, covariates=~ollwsdrt, subject=ptnum, data=psor, qmatrix = psor.q)
get.covhist(psor1.msm)
observed.msm(psor1.msm)
expected.msm(psor1.msm)
expected.msm(psor1.msm, covariates="mean")
}, NA)
})
test_that("observed, expected etc with PCI, should run",{
expect_error({
psor1.msm <- msm(state ~ months, covariates=~ollwsdrt+hieffusn, subject=ptnum, data=psor, qmatrix = psor.q, pci=c(5,10))
covhist <- get.covhist(psor1.msm)
observed.msm(psor1.msm)
expected.msm(psor1.msm)
expected.msm(psor1.msm, covariates="mean")
}, NA)
})
test_that("subset argument to observed",{
subs <- psor.msm$data$mf$"(subject)"[!duplicated(psor.msm$data$mf$"(subject)") & psor.msm$data$mf$ollwsdrt==0]
expect_equivalent(observed.msm(psor.msm, subset=subs)$obstab[,"State 4"], c(0, 5, 20, 30, 37, 38, 40, 41, 43, 44, 45))
})
test_that("error handling: formula",{
## formula
expect_error(msm(), "state ~ time formula not given")
expect_error(cav.msm <- msm(state, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "not found")
expect_error(cav.msm <- msm(~1, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "invalid data")
expect_error(cav.msm <- msm("foo", subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "not a formula")
})
test_that("error handling: qmatrix",{
wrong.q <- cbind(c(0,1,2), c(0,1,2))
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"Number of rows and columns of qmatrix should be equal")
wrong.q <- cbind(c(0,1), c(0,1))
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, fixedpars=TRUE),"State vector contains elements not in 1, 2")
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"Not all the states specified in \"deathexact\" are absorbing")
wrong.q <- "foo"
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"qmatrix should be a numeric matrix")
wrong.q <- 1
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"qmatrix should be a numeric matrix")
})
test_that("error handling: subject",{
expect_error(cav.msm <- msm(state~years, subject="foo", data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),
"variable lengths differ")
expect_error(cav.msm <- msm(state~years, subject=foo, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),
"not found")
})
test_that("error handling: obstype",{
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype="foo", deathexact = TRUE, fixedpars=TRUE),"should be numeric")
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=rep(1,10), deathexact = TRUE, fixedpars=TRUE),"obstype of length") ##FIXME
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=rep(4, nrow(cav)), deathexact = TRUE, fixedpars=TRUE),"elements of obstype should be 1, 2, or 3")
obstype <- rep(1, nrow(cav))
obstype[c(1,8)] <- 5
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=obstype, deathexact = TRUE, fixedpars=TRUE), NA) # no error: obstype for first subject doesn't matter
obstype[2] <- 5
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=obstype, deathexact = TRUE, fixedpars=TRUE),"elements of obstype should be 1, 2, or 3") # error
})
test_that("error handling: covariates",{
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = "wibble"),"should be a formula or list of formulae")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sux),"not found")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=TRUE, misccovariates = "wobble"),"should be a formula")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sox),"not found")
})
test_that("error handling: covinits",{
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits="foo", fixedpars=TRUE),"covinits should be a list")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(sex="foo", age="bar"), fixedpars=TRUE),"should be numeric")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(sex=c(1,2,3), age="bar"), fixedpars=TRUE),"should be a list of numeric vectors")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(age=rep(0.1, 7)), fixedpars=TRUE),"covariate age in covinits unknown")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(age=rep(0.1, 7), foo=1, bar=2), fixedpars=TRUE),"covariates age, foo, bar in covinits unknown")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(sex=1), fixedpars=TRUE),"initial values of length 1, should be 7")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, misccovinits="foo", fixedpars=TRUE), "hcovinits should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, misccovariates = ~ sex, misccovinits=list(sex=1, age="bar"), fixedpars=TRUE), "misccovariates supplied but no ematrix")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, misccovariates = ~ sex, misccovinits=list(sex=1, age="bar"), fixedpars=TRUE),"misccovariates supplied but no ematrix")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, ematrix=ematrix, misccovariates = ~ sex, misccovinits=list(sex=1, age="bar"), fixedpars=TRUE),"covinits should be a list of numeric")
})
test_that("error handling: constraints",{
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, constraint="foo"),"constraint should be a list")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, constraint=list(foo="bar")),"constraint should be a list of numeric vectors")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, constraint=list(foo=1)),"Covariate .+ in constraint statement not in model.")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, constraint="foo", fixedpars=TRUE),"constraint specified but no covariates")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, miscconstraint=list(foo="bar")),"constraint should be a list of numeric vectors")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, miscconstraint=list(foo=1)),"Covariate .+ in constraint statement not in model.")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, miscconstraint=list(sex=1)),"constraint of length 1, should be 4")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, qconstraint="foo", fixedpars=TRUE),"qconstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, qconstraint=list(c(1,1,2)), fixedpars=TRUE),"qconstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, qconstraint=c(1,1,2), fixedpars=TRUE),"baseline intensity constraint of length 3, should be 7")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, econstraint="foo", fixedpars=TRUE) ,"econstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, econstraint=list(c(1,1,2)), fixedpars=TRUE) ,"econstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, econstraint=c(1,1,2), fixedpars=TRUE),"baseline misclassification constraint of length 3, should be 4")
})
test_that("error handling: initprobs",{
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, initprobs="poo", fixedpars=TRUE),"initprobs should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, initprobs=c(1,2), fixedpars=TRUE),"initprobs vector of length 2, should be vector of length 4 or a matrix")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, initprobs=c(2,1,1,1), fixedpars=TRUE), NA) # scaled to sum to 1.
})
test_that("error handling: check states",{
wrong.q <- cbind(c(0,1), c(0,1)) # extra states in data
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, fixedpars=TRUE),"State vector contains elements not in 1, 2")
wrong.q <- rbind(c(0,1,2,3,1), c(0,1,3,4,1), c(0,1,2,3,2), c(0,1,2,3,4), c(0,0,0,0,0))
expect_warning(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, fixedpars=TRUE),"State vector doesn't contain observations of 5")
})
test_that("error handling: check times",{
cav.wrong <- cav
cav.wrong$years[3:5] <- 4:2
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),"not ordered by time")
cav.wrong <- cav
cav.wrong$PTNUM[4:5] <- 100003
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),"Observations within subjects .+ are not adjacent in the data")
cav2 <- cav
cav2$years[10] <- cav2$years[9]
expect_warning(msm(state ~ years, subject=PTNUM, data = cav2, qmatrix = twoway4.q, fixedpars=TRUE),
"Different states observed at the same time on the same subject at observations 9 and 10")
## with missing data
cav2$years[6] <- NA ## report original rows before excluding missing data
expect_warning(msm(state ~ years, subject=PTNUM, data = cav2, qmatrix = twoway4.q, fixedpars=TRUE),
"Different states observed at the same time on the same subject at observations 9 and 10")
cav2 <- cav2[6:nrow(cav),]
expect_warning(msm(state ~ years, subject=PTNUM, data = cav2, qmatrix = twoway4.q, fixedpars=TRUE),
"Different states observed at the same time on the same subject at observations 4 and 5",
"Subject 100002 only has one complete observation")
})
test_that("error handling: check model",{
cav.wrong <- cav
cav.wrong$state[4] <- 1
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = oneway4.q, deathexact = TRUE, fixedpars=TRUE),"Data may be inconsistent with transition matrix for model without misclassification:\nindividual 100002 moves from state 2 to state 1 at observation 4")
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.i, gen.inits=TRUE, fixedpars=TRUE),"individual 100046 moves from state 2 to state 1 at observation 225")
## row number reporting with missing data - should be the same
cav.wrong$PTNUM[2] <- NA
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = oneway4.q, deathexact = TRUE, fixedpars=TRUE),"100002 moves from state 2 to state 1 at observation 4")
cav.wrong <- cav
cav.wrong$state[4] <- 1
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, exacttimes=TRUE, fixedpars=TRUE),"individual 100003 moves from state 1 to state 3 at observation 10")
## row number reporting with missing data
cav.wrong$state[3] <- NA
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, exacttimes=TRUE, fixedpars=TRUE),"individual 100003 moves from state 1 to state 3 at observation 10") # should complain about obs 10
obstype <- rep(2, nrow(cav))
obstype[10] <- 1
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, obstype=obstype, fixedpars=TRUE),"individual 100006 moves from state 1 to state 3 at observation 29")
## absorbing-absorbing transitions
cav.wrong$state[6] <- 4
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, obstype=obstype, fixedpars=TRUE),"Absorbing - absorbing transition at observation 7")
})
test_that("error handling: death",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = "foo", fixedpars=TRUE) ,"Exact death states indicator must be numeric")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = 5, fixedpars=TRUE) ,"Exact death states indicator contains states not in 1, 2, ... , 4")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = 1:5, fixedpars=TRUE) ,"Exact death states indicator contains states not in 1, 2, ... , 4")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = 3, fixedpars=TRUE),"Not all the states specified in \"deathexact\" are absorbing" )
})
test_that("error handling: censor",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, censor="rubbish", fixedpars=TRUE),"censor must be numeric")
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, censor=1, fixedpars=TRUE),"some censoring indicators are the same as actual states")
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, censor.states="rubbish", fixedpars=TRUE) ,"censor.states supplied but censor not supplied")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.cens, qmatrix = twoway4.q, censor=99, censor.states="rubbish", fixedpars=TRUE) ,"censor.states should be all numeric")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.cens, qmatrix = twoway4.q, censor=99, censor.states=list(c(1,2,3), "rubbish"), fixedpars=TRUE) ,"censor.states should be a vector")
})
test_that("error handling: obstype",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype="rubbish", fixedpars=TRUE) ,"obstype should be numeric")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=c(1,2,3), fixedpars=TRUE) ,"obstype of length 3, should be length 1 or 2846")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=4, fixedpars=TRUE),"elements of obstype should be 1, 2, or 3" )
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=rep(4,nrow(cav)), fixedpars=TRUE) ,"elements of obstype should be 1, 2, or 3")
})
test_that("error handling: fixedpars",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars="foo"),"Elements of fixedpars should be in 1, ..., 7")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=list(c(1,3,4))),"Elements of fixedpars should be in 1, ..., 7")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=1:8),"Elements of fixedpars should be in 1, ..., 7")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=0.5),"Elements of fixedpars should be in 1, ..., 7")
})
test_that("error handling: plot",{
expect_error(plot.msm("foo"),"expected .+ msm model")
expect_error(plot.msm(psor.msm, from="foo"), "from must be numeric")
expect_error(plot.msm(psor.msm, to="foo"), "to must be numeric")
expect_error(plot.msm(psor.msm, from = 1:8, to=3),"from must be a vector of states in 1, ..., 4")
expect_error(plot.msm(psor.msm, to = 3),"to must be an absorbing state")
expect_error(plot.msm(psor.msm, range="foo"))
expect_error(plot.msm(psor.msm, range=1:6),"range must be a numeric vector of two elements")
})
test_that("recreate.olddata",{
expect_error(recreate.olddata(psor.msm), NA)
})
test_that("form.output",{
qo <- msm.form.qoutput(psor.msm)
expect_equal(qo$base.Estimate[2],
qmatrix.msm(psor.msm, covariates="mean")$estimates[1,2])
print(qmatrix.msm(psor.msm))
})
test_that("plots",{
skip_on_cran()
expect_error({
plot.msm(psor.msm)
plotprog.msm(state ~ months, subject=ptnum, data=psor)
plot.survfit.msm(psor.msm)
contour.msm(psor.msm)
persp.msm(psor.msm)
image.msm(psor.msm)
surface.msm(psor.msm, xrange=c(-2.4,-2.3))
}, NA)
})
test_that("miscellaneous",{
expect_error(msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~ollwsdrt+hieffusn,
constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), control=list(fnscale=1),
na.action=na.fail))
expect_error(model.matrix(psor.msm), NA)
expect_error(model.frame(psor.msm), NA)
})
test_that("single-column matrix in covariates",{
psor.q <- rbind(c(0,0.1,0,0),c(0,0,0.1,0),c(0,0,0,0.1),c(0,0,0,0))
psor$matcov <- matrix(psor$ollwsdrt, ncol=1)
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~matcov, fixedpars=TRUE)
psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~ollwsdrt, fixedpars=TRUE)
expect_equal(psor.msm$minus2loglik, psor2.msm$minus2loglik)
})
Try the msm package in your browser
Any scripts or data that you put into this service are public.
msm documentation built on Nov. 24, 2023, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
context("msm simple model likelihoods")
test_that("simple model, exact death times",{
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
method="BFGS", control=list(trace=5, REPORT=1))
print(cav.msm)
printold.msm(cav.msm)
expect_equal(4908.81676837903, cav.msm$minus2loglik)
})
test_that("simple model, not exact death times",{
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = FALSE, fixedpars=TRUE)
expect_equal(4833.00640639644, cav.msm$minus2loglik)
})
test_that("simple model, covariates",{
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
covariates = ~ sex, covinits = list(sex=rep(0.01, 7)))
expect_equal(4909.17147259115, cav.msm$minus2loglik)
})
test_that("autogenerated inits reproduce crudeinits",{
cinits <- crudeinits.msm(state ~ years, PTNUM, data=cav, qmatrix=twoway4.q)
expect_equal(msm( state ~ years, subject=PTNUM, data = cav, qmatrix = cinits, deathexact = TRUE, fixedpars=TRUE)$minus2loglik,
msm( state ~ years, subject=PTNUM, data = cav, qmatrix=twoway4.i, gen.inits=TRUE, deathexact = TRUE, fixedpars=TRUE)$minus2loglik)
})
test_that("data as global variables",{
state.g <- cav$state; time.g <- cav$years; subj.g <- cav$PTNUM
cav.msm <- msm(state.g ~ time.g, subject=subj.g, qmatrix = twoway4.i, gen.inits=TRUE, fixedpars=TRUE)
expect_equal(4119.9736299032, cav.msm$minus2loglik)
})
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), control=list(fnscale=1))
test_that("psor model: covariates, constraints",{
print(psor.msm)
printold.msm(psor.msm)
expect_equal(1114.89946121717, psor.msm$minus2loglik, tol=1e-06)
expect_equal(0.0959350004999946, psor.msm$Qmatrices$baseline[1,2], tol=1e-06)
expect_equal(exp(psor.msm$Qmatrices$logbaseline[c(5,10,15)]), psor.msm$Qmatrices$baseline[c(5,10,15)], tol=1e-06)
})
test_that("psor model: transition-specific covariates",{
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, covariates = list("2-3"=~ollwsdrt, "3-4"=~hieffusn), control=list(fnscale=1))
expect_equal(hazard.msm(psor.msm)$ollwsdrt["State 1 - State 2","HR"], 1)
expect_equal(hazard.msm(psor.msm)$hieffusn["State 2 - State 3","HR"], 1)
expect_error(msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = list("meep"=~ollwsdrt, "3-4"=~hieffusn), control=list(fnscale=1)),
"not in format \"number-number\"")
})
test_that("psor model: transition-specific covariates with pci",{
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = list("2-3"=~ollwsdrt, "3-4"=~hieffusn),
pci = 10, fixedpars=7:8, control=list(fnscale=1))
h <- hazard.msm(psor.msm)
expect_equal(h$ollwsdrt["State 1 - State 2","HR"], 1)
expect_equal(h$hieffusn["State 1 - State 2","HR"], 1)
expect_true(!isTRUE(all.equal(h$"timeperiod[10,Inf)"["State 1 - State 2","HR"], 1)))
expect_equal(h$"timeperiod[10,Inf)"["State 2 - State 3","HR"], 1)
expect_equal(h$"timeperiod[10,Inf)"["State 3 - State 4","HR"], 1)
})
psor.nocen.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), center=FALSE)
test_that("no covariate centering",{
expect_equal(exp(psor.nocen.msm$Qmatrices$logbaseline[c(5,10,15)]), psor.nocen.msm$Qmatrices$baseline[c(5,10,15)], tol=1e-06)
})
test_that("gen.inits with missing values for state / time",{
psor2 <- psor; psor2$ptnum[13:14] <- psor2$months[7:8] <- psor2$state[7:8] <- NA
psor2.msm <- msm(state ~ months, subject=ptnum, data=psor2, gen.inits=TRUE, qmatrix = psor.q, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), fixedpars=TRUE)
expect_equal(1179.95441284169, psor2.msm$minus2loglik, tol=1e-06)
})
test_that("exact transition times using exacttimes",{
msmtest5 <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, exacttimes=TRUE, fixedpars=TRUE)
expect_equal(3057.85781916437, msmtest5$minus2loglik, tol=1e-06)
})
test_that("exact transition times using obstype vector of all 2",{
msmtest5 <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, obstype=rep(2, nrow(bos)), fixedpars=TRUE)
expect_equal(3057.85781916437, msmtest5$minus2loglik, tol=1e-06)
})
test_that("exact transition times with death, should be same",{
expect_warning(msmtested <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, deathexact=5, obstype=rep(2, nrow(bos)), exacttimes=TRUE, fixedpars=TRUE), "Ignoring death")
expect_equal(3057.85781916437, msmtested$minus2loglik, tol=1e-06)
(msmtest5 <- msm(state ~ time, qmatrix = fiveq, subject = ptnum, data = bos, deathexact=5, obstype=rep(2, nrow(bos)), fixedpars=TRUE)) # no warning, inconsistently
expect_equal(3057.85781916437, msmtest5$minus2loglik, tol=1e-06)
})
cav$statefac <- factor(cav$state)
cav$statefac2 <- factor(cav$state, labels=c("none","mild","severe","death"))
test_that("factors as states, death",{
cav.msm <- msm( statefac ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
method="BFGS", control=list(trace=5, REPORT=1))
expect_equal(4908.81676837903, cav.msm$minus2loglik)
expect_error(msm( statefac2 ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "state variable should be numeric or a factor with ordinal numbers as levels")
})
test_that("factor covariates with factor() in formula",{
## Should be no need for users to do this. factors should already be identified as such in the data frame
expect_warning(cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ factor(pdiag), covinits=list(sex=rep(0.1,7)), fixedpars=TRUE), "covariate .+ unknown")
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ factor(pdiag), covinits=list("factor(pdiag)Hyper"=rep(0.1,7)), fixedpars=TRUE)
expect_equal(4793.20440566637, cavfaccov.msm$minus2loglik, tol=1e-06)
})
test_that("factors as global variables",{
state.g <- cav$state; time.g <- cav$years; subj.g <- cav$PTNUM; pdiag.g <- factor(cav$pdiag)
cavfaccov.msm <- msm(state.g ~ time.g, subject=subj.g, qmatrix = twoway4.q, covariates = ~ pdiag.g, fixedpars=TRUE)
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
})
test_that("factor covariates using existing factors: inits are given to contrasts ",{
## Warnings could be more informative here
expect_warning(cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ pdiag, covinits=list(pdiag=rep(0.1,7)), fixedpars=TRUE), "covariate .+ unknown")
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
expect_warning(cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ pdiag, covinits=list(pdiagNonexistentlevel=rep(0.1,7)), fixedpars=TRUE), "covariate .+ unknown")
expect_equal(4793.30238295565, cavfaccov.msm$minus2loglik, tol=1e-06)
cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q,
covariates = ~ pdiag, covinits=list(pdiagHyper=rep(0.1,7)), fixedpars=TRUE)
expect_equal(4793.20440566637, cavfaccov.msm$minus2loglik, tol=1e-06)
cavfaccov.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ pdiag,
covinits=list(pdiagHyper=rep(0.1,7),pdiagIDC=rep(0.1,7),pdiagIHD=rep(0.1,7),pdiagOther=rep(0.1,7),pdiagRestr=rep(0.1,7)), fixedpars=TRUE) # OK
expect_equal(4793.13035886505, cavfaccov.msm$minus2loglik, tol=1e-06)
})
context("censored states")
test_that("censored states: final state censored",{
cavcens.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens, qmatrix=twoway4.q, censor=99, fixedpars=TRUE)
expect_equal(4724.26606344485, cavcens.msm$minus2loglik, tol=1e-06)
v <- viterbi.msm(cavcens.msm)
expect_equal(v$observed[v$observed<10], v$fitted[v$observed<10])
expect_equal(v$fitted[v$observed==99][1], 3)
})
test_that("two kinds of censoring",{
cavcens2.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens2, qmatrix=twoway4.q, censor=c(99, 999), censor.states=list(c(1,2,3), c(2,3)), fixedpars=TRUE)
expect_equal(4678.23348518727, cavcens2.msm$minus2loglik, tol=1e-06)
v <- viterbi.msm(cavcens2.msm)
expect_equal(v$observed[v$observed<10], v$fitted[v$observed<10])
expect_equal(v$fitted[v$observed==99][1], 3)
})
test_that("intermediate state censored",{
cavcens3.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens3, qmatrix=twoway4.q, censor=c(99, 999), censor.states=list(c(2,3), c(1,2,3)), fixedpars=TRUE)
expect_equal(4680.66073438518, cavcens3.msm$minus2loglik, tol=1e-06)
v <- viterbi.msm(cavcens3.msm)
expect_equal(v$observed[v$observed<10], v$fitted[v$observed<10])
expect_true(all(v$fitted[v$observed==99] %in% 2:3))
expect_true(all(v$fitted[v$observed==999] %in% 1:3))
})
test_that("first state censored",{
cav.cens4 <- cav
cav.cens4$state[c(1,8,12,22)] <- 99
cavcens4.msm <- msm(state ~ years, subject=PTNUM, data=cav.cens4, qmatrix=twoway4.q,
censor=c(99), censor.states=list(c(2,3)), fixedpars=TRUE)
expect_equal(4846.06045097812, cavcens4.msm$minus2loglik)
v <- viterbi.msm(cavcens4.msm)
expect_true(all(v$fitted[v$observed==99] %in% 2:3))
})
test_that("piecewise constant intensities with pci",{
cav5.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE,
pci = c(5), covinits = list("timeperiod[5,Inf)"=rep(0.01,7)),
)
expect_equal(4906.01423796805, cav5.msm$minus2loglik, tol=1e-06)
cav10.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, pci = c(5,10), fixedpars=TRUE,
covinits = list("timeperiod[5,10)"=rep(0.01,7), "timeperiod[10,Inf)"=rep(0.01,7)),
)
expect_equal(4905.61646158639, cav10.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with pci, cut points outside data",{
## Make sure works for pci outside time range, with warning
expect_warning(cav5.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE, pci = c(-1,5,50,60), covinits = list("timeperiod[5,Inf)"=rep(0.01,7))), "cut point .+ less than")
expect_warning(cav.pci.msm <- msm( state ~ years, subject=PTNUM, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE, pci = c(-1,50,60)))
cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE)
expect_equal(cav.pci.msm$minus2loglik, cav.msm$minus2loglik) # degrades to time-homogeneous if all cuts outside data
})
suppressWarnings(RNGversion("3.5.0"))
set.seed(22061976)
cav$pdiag3 <- cav$pdiag
cav$pdiag3[!cav$pdiag %in% c("IDC","IHD")] <- "Other"
cav$pdiag3 <- factor(cav$pdiag3)
subs <- cav$PTNUM %in% sample(unique(cav$PTNUM), 50)
cavsub <- subset(cav, subs)
cavsub$maxtime <- tapply(cavsub$years, cavsub$PTNUM, max)[as.character(cavsub$PTNUM)]
cavsub.extra <- cavsub[cavsub$years==0 & cavsub$maxtime >= 5,]
cavsub.extra$years <- 5
cavsub.extra$state <- 99
cavsub.extra10 <- cavsub[cavsub$years==0 & cavsub$maxtime >= 10,]
cavsub.extra10$years <- 10
cavsub.extra10$state <- 999
cavsub2 <- rbind(cavsub, cavsub.extra, cavsub.extra10)
cavsub2 <- cavsub2[order(cavsub2$PTNUM, cavsub2$years),]
cavsub2$after5 <- ifelse(cavsub2$years>=5 & cavsub2$years<10, 1, 0)
cavsub2$after10 <- ifelse(cavsub2$years>=10, 1, 0)
cavsub2$after510 <- as.numeric(cavsub2$after5 | cavsub2$after10)
test_that("piecewise constant intensities with pci, with other covariates",{
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub,
qmatrix = twoway4.q, covariates = ~ pdiag3 + sex, deathexact = TRUE, pci = c(5,10), fixedpars=TRUE,
covinits = list("timeperiod[5,10)"=rep(0.01,7), "timeperiod[10,Inf)"=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
)
expect_equal(448.122802051545, cav5cov.msm$minus2loglik, tol=1e-06)
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub2,
qmatrix = twoway4.q, deathexact = TRUE, covariates = ~ pdiag3 + sex + after5 + after10,
censor=c(99,999), censor.states=list(1:4,1:4),
covinits = list(after5=rep(0.01,7), after10=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
fixedpars=TRUE
)
expect_equal(448.122802051545, cav5cov.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with pci, with uncentered covariates",{
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub,
qmatrix = twoway4.q, covariates = ~ pdiag3 + sex, deathexact = TRUE, pci = c(5,10), fixedpars=TRUE, center=FALSE,
covinits = list("timeperiod[5,10)"=rep(0.01,7), "timeperiod[10,Inf)"=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
)
expect_equal(449.691983558378, cav5cov.msm$minus2loglik, tol=1e-06)
cav5cov.msm <- msm( state ~ years, subject=PTNUM, data = cavsub2,
qmatrix = twoway4.q, deathexact = TRUE, covariates = ~ pdiag3 + sex + after5 + after10,
censor=c(99,999), censor.states=list(1:4,1:4), center=FALSE,
covinits = list(after5=rep(0.01,7), after10=rep(0.01,7), pdiag3IHD=rep(0.01,7), pdiag3Other=rep(0.01,7)),
fixedpars=TRUE
)
expect_equal(449.691983558378, cav5cov.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with pci, with other censored states",{
cav.cens <- cav
cav.cens$state[cav.cens$state==4][1:50] <- 99
cav5cens.msm <- msm(state ~ years, subject=PTNUM, data = cav.cens,
qmatrix = twoway4.q, deathexact = TRUE, censor=99, pci = 5,
covinits=list("timeperiod[5,Inf)"=rep(0.02,7)), fixedpars=TRUE,
method="BFGS", control=list(trace=5, REPORT=1))
expect_equal(4754.41981265159, cav5cens.msm$minus2loglik, tol=1e-06)
})
test_that("piecewise constant intensities with covariates in HMMs",{
misccov.msm <- msm(state ~ years, subject = PTNUM, data = cav,
qmatrix = oneway4.q, ematrix=ematrix, deathexact = 4, fixedpars=TRUE,
pci = 5, covinits=list("timeperiod[5,Inf)"=rep(0.0001,5)),
misccovariates = ~dage + sex, misccovinits = list(dage=c(0.01,0.02,0.03,0.04), sex=c(-0.013,-0.014,-0.015,-0.016)))
expect_equal(4306.29865288466, misccov.msm$minus2loglik, tol=1e-06)
})
context("output functions")
test_that("qmatrix.msm for psor model, defaults",{
expect_equal(c(-0.0959350004999946, 0, 0, 0, 0.0959350004999946, -0.164306508892574, 0, 0, 0, 0.164306508892574, -0.254382807485639, 0, 0, 0, 0.254382807485639, 0), as.numeric(qmatrix.msm(psor.msm)$estimates), tol=1e-03)
expect_equal(c(0.0115942726096754, 0, 0, 0, 0.0115942726096754, 0.0196169975000406, 0, 0, 0, 0.0196169975000406, 0.0375066077515386, 0, 0, 0, 0.0375066077515386, 0), as.numeric(qmatrix.msm(psor.msm)$SE), tol=1e-03)
expect_error(qmatrix.msm(psor.msm, ci="normal", cl=-1), "expected cl")
expect_equivalent(qmatrix.msm(psor.msm, sojourn=TRUE)$sojourn[1:3],
sojourn.msm(psor.msm)$estimates)
qmatrix.msm(psor.msm, ci="normal", B=2)
qmatrix.msm(psor.msm, sojourn=TRUE, ci="normal", B=2)
expect_error(qmatrix.msm("foo"), "expected .+ msm model")
expect_error(qmatrix.msm(psor.msm, covariates="foo"), "covariates argument must be")
expect_warning(qmatrix.msm(psor.msm, covariates=list(foo=1)), "ignoring")
})
test_that("qmatrix.msm defaults to Qmatrices in object",{
expect_equal(psor.msm$Qmatrices$baseline, qmatrix.msm(psor.msm, covariates="mean", ci="none"))
expect_equal(psor.nocen.msm$Qmatrices$baseline, qmatrix.msm(psor.nocen.msm, covariates=0, ci="none"))
})
test_that("qmatrix.msm with supplied covariates",{
qmat <- qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4))
expect_equal(c(-0.121430585652200, 0, 0, 0, 0.121430585652200, -0.207972362475868, 0, 0, 0, 0.207972362475868, -0.257535341208494, 0, 0, 0, 0.257535341208494, 0), as.numeric(qmat$estimates), tol=1e-03)
expect_equal(c(0.0162156605802465, 0, 0, 0, 0.0162156605802465, 0.0266727053124233, 0, 0, 0, 0.0266727053124233, 0.0364321127089265, 0, 0, 0, 0.0364321127089265, 0), as.numeric(qmat$SE), tol=1e-04)
expect_warning(qmatrix.msm(psor.msm, covariates=list(hieffusn=0.1, foo=0.4)), "Covariate .+ unknown")
})
test_that("qmatrix.msm with non-default confidence limits",{
qmat <- qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4), cl=0.99)
expect_equal(c(-0.171282667596986, 0, 0, 0, 0.0860880282792585, -0.289385121267802, 0, 0, 0, 0.149463467106753, -0.370756460718086, 0, 0, 0, 0.178889538008097, 0), as.numeric(qmat$L), tol=1e-04)
})
test_that("qmatrix.msm sojourn component",{
soj <- qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4), sojourn=TRUE)$sojourn
expect_equal(c(8.23515751512713, 4.80833120370037, 3.88296221911705, Inf), as.numeric(soj), tol=1e-03)
expect_equal(as.numeric(soj[1:3]), sojourn.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4))[,"estimates"])
})
test_that("qmatrix.msm bug for user-supplied covariates fixed in 1.1.3",{
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$SE[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0, ollwsdrt=0))$SE[1,2])
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$SE[1,1],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0, ollwsdrt=0))$SE[1,1])
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$L[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0, ollwsdrt=0))$L[1,2])
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=0)$L[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0))$L[1,2])
cm <- psor.nocen.msm$qcmodel$covmeans
expect_equal(qmatrix.msm(psor.nocen.msm, covariates="mean")$SE[1,2],
qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=cm["hieffusn"], ollwsdrt=cm["ollwsdrt"]))$SE[1,2])
})
test_that("qmatrix.msm: unspecified covariate values default to zero",{
expect_equal(psor.nocen.msm$Qmatrices$baseline, qmatrix.msm(psor.nocen.msm, covariates=list(ollwsdrt=0), ci="none")) # missing covs default to zero
expect_equal(qmatrix.msm(psor.msm, covariates=list(hieffusn=0)),
qmatrix.msm(psor.msm, covariates=list(ollwsdrt=0, hieffusn=0)))
expect_equal(qmatrix.msm(psor.nocen.msm, covariates=list(hieffusn=0)),
qmatrix.msm(psor.nocen.msm, covariates=list(ollwsdrt=0, hieffusn=0)))
})
test_that("sojourn.msm",{
expect_equal(psor.msm$sojourn, sojourn.msm(psor.msm, covariates="mean"))
expect_equal(psor.nocen.msm$sojourn, sojourn.msm(psor.nocen.msm, covariates=0))
soj <- sojourn.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4))
expect_equal(c(8.23515751512713, 4.80833120370037, 3.88296221911705, 1.09971073904434, 0.616674252838334, 0.549301375677405, 6.33875136203292, 3.73961380505919, 2.94271599303942, 10.6989240349703, 6.18246967994404, 5.12363260020806), as.numeric(unlist(soj)), tol=1e-04)
soj <- sojourn.msm(psor.msm, covariates=list(ollwsdrt=0.1, hieffusn=0.4), cl=0.99)
expect_equal(5.83830234564607, soj[1,"L"], tol=1e-04)
expect_error(sojourn.msm("foo"), "expected .+ msm model")
expect_error(sojourn.msm(psor.msm, covariates="foo"),"covariates argument must be")
expect_warning(sojourn.msm(psor.msm, covariates=list(foo=1)), "ignoring")
})
test_that("pmatrix.msm",{
expect_equal(0.149287738928777, pmatrix.msm(psor.msm, ci="none", t=10)[1,3], tol=1e-04)
p <- pmatrix.msm(psor.msm, t=10, covariates=list(ollwsdrt=0.1, hieffusn=0.2))
expect_equal(0.18196160265907, p[1,3], tol=1e-04)
set.seed(22061976); expect_equal(0.12, pmatrix.msm(psor.msm, ci="normal", B=3)$L[2,3], tol=1e-01, scale=1)
expect_error(pmatrix.msm("foo"), "expected .+ msm model")
expect_error(pmatrix.msm(psor.msm, t="foo"), "must be a positive number")
expect_error(pmatrix.msm(psor.msm, -9), "must be a positive number")
expect_equivalent(unclass(pmatrix.msm(psor.msm, 0)), diag(4))
expect_warning(pmatrix.msm(psor.msm, 1, covariates=list(foo=1)), "ignoring")
})
test_that("pnext.msm",{
expect_equal(pnext.msm(psor.msm, ci="none")$estimate[1,2], 1)
expect_equal(pnext.msm(psor.msm, ci="delta")$estimate[1,2], 1)
expect_equal(pnext.msm(psor.msm, ci="normal", B=2)$estimate[1,2], 1)
})
test_that("qratio.msm",{
q <- qratio.msm(psor.msm, c(1,2), c(2,3))
expect_equal(c(0.583878474075081, 0.0996029045389022, 0.417943274168735, 0.815694601537263), as.numeric(q), tol=1e-04)
q <- qratio.msm(psor.msm, c(1,2), c(2,3), cl=0.99)
expect_equal(0.376262194364283, as.numeric(q["L"]), tol=1e-04)
q <- qratio.msm(psor.msm, c(1,1), c(2,3))
expect_equal(c(-0.583878474075081, 0.0996029045389022, -0.815694601537263, -0.417943274168735), as.numeric(q), tol=1e-04)
q <- qratio.msm(psor.msm, c(2,2), c(2,3))
expect_equivalent(c(-1,0,-1,-1), q)
qratio.msm(psor.msm, c(1,2), c(2,3), ci="norm", B=2)
expect_error(qratio.msm("foo"))
expect_error(qratio.msm(psor.msm, "foo"))
expect_error(qratio.msm(psor.msm, c(1,8), c(1,0)))
expect_error(qratio.msm(psor.msm, c(1,2), c(1,0)))
expect_error(qratio.msm(psor.msm, c(1,2), c(2,3), cl="foo"))
expect_error(qratio.msm(psor.msm, c(1,2), c(2,3), cl=2), "expected cl in")
qq <- qratio.msm(psor.msm, c(1,2), c(2,2))
QQ <- qmatrix.msm(psor.msm)
expect_equal(qq[["estimate"]], (QQ[1,2]/QQ[2,2])[["estimate"]])
qq <- qratio.msm(psor.msm, c(1,1), c(2,2))
QQ <- qmatrix.msm(psor.msm)
expect_equal(qq[["estimate"]], (QQ[1,1]/QQ[2,2])[["estimate"]])
})
test_that("coef.msm",{
co <- coef.msm(psor.msm)
expect_equal(0.498319866154661, co$hieffusn[1,2], tol=1e-04)
expect_error(coef.msm("foo"))
})
test_that("hazard.msm",{
haz <- hazard.msm(psor.msm)
expect_equal(0.385347226135311, haz$ollwsdrt[1,2], tol=1e-04)
expect_equal(0.385347226135311, haz$ollwsdrt[2,2], tol=1e-04)
expect_equal(2.35928404626333, haz$hieffusn[1,3], tol=1e-04)
expect_equal(2.35928404626333, haz$hieffusn[3,3], tol=1e-04)
haz <- hazard.msm(psor.msm, hazard.scale=2)
expect_equal(0.148492484690178, haz$ollwsdrt[1,2], tol=1e-04)
expect_equal(haz$ollwsdrt[1,2], haz$ollwsdrt[2,2])
expect_equal(haz$hieffusn[1,3], haz$hieffusn[3,3])
haz <- hazard.msm(psor.msm, hazard.scale=c(1,2))
expect_equal(0.385347226135311, haz$ollwsdrt[1,2], tol=1e-04)
expect_equal(haz$ollwsdrt[1,2], haz$ollwsdrt[2,2], tol=1e-04)
expect_equal(haz$hieffusn[1,3], haz$hieffusn[3,3])
expect_error(hazard.msm("foo"))
expect_error(hazard.msm(psor.msm, hazard.scale="foo"))
expect_error(hazard.msm(psor.msm, hazard.scale=c(1,2,3)), "hazard.scale of length")
})
test_that("transient.msm and absorbing.msm",{
expect_equivalent(c(1,2,3), transient.msm(psor.msm))
expect_equivalent(4, absorbing.msm(psor.msm))
expect_error(transient.msm("foo"))
expect_error(absorbing.msm("foo"))
expect_error(transient.msm(qmatrix="foo"))
expect_error(absorbing.msm(qmatrix="foo"))
expect_error(transient.msm(qmatrix=c(1,4,5,6)))
expect_error(absorbing.msm(qmatrix=c(1,4,5,6)))
expect_error(transient.msm(qmatrix=cbind(c(1,2,3),c(1,3,2))))
expect_error(absorbing.msm(qmatrix=cbind(c(1,2,3),c(1,3,2))))
expect_error(transient.msm())
expect_error(absorbing.msm())
})
test_that("prevalence.msm",{
p <- prevalence.msm(psor.msm)
expect_equal(59, p$Observed[5,5], tol=1e-06)
expect_equal(59, p$Expected[5,5], tol=1e-06)
expect_equal(57.35294, p$"Observed percentages"[4,4], tol=1e-03)
expect_equal(49.96882, p$"Expected percentages"[4,4], tol=1e-03)
summ <- summary.msm(psor.msm)
expect_equal(summ$prevalences, p)
p <- prevalence.msm(psor.msm, times=seq(0,60,5))
expect_equal(63, p$Observed[5,5], tol=1e-06)
expect_equal(63, p$Expected[5,5], tol=1e-06)
expect_equal(50.70423, p$"Observed percentages"[4,4], tol=1e-03)
expect_equal(41.46338, p$"Expected percentages"[4,4], tol=1e-03)
expect_error(prevalence.msm("foo"))
expect_error(summary.msm("foo"))
p <- prevalence.msm(psor.msm, covariates=list(hieffusn=0, ollwsdrt=1))
expect_equal(p$Observed[1,1], 1)
p <- prevalence.msm(psor.msm, covariates=list(hieffusn=0, ollwsdrt=1),
ci="normal", B=10)
expect_true(is.numeric(p[[2]]$ci[,,"2.5%"][1,1]))
plot.prevalence.msm(psor.msm)
})
test_that("pmatrix.piecewise.msm",{
times <- c(5, 10, 15)
covariates <- list(list(ollwsdrt=0, hieffusn=0),
list(ollwsdrt=0, hieffusn=1),
list(ollwsdrt=1, hieffusn=0),
list(ollwsdrt=1, hieffusn=1)
)
p <- pmatrix.msm(psor.msm, 3, covariates=covariates[[1]])
pp <- pmatrix.piecewise.msm(psor.msm, 0, 3, times, covariates)
expect_equal(pp[1,3], p[1,3], tol=1e-04)
p <- pmatrix.piecewise.msm(psor.msm, 0, 7, times, covariates)
expect_equal(0.172773087945103, p[1,3], tol=1e-04)
pp <- pmatrix.piecewise.msm(psor.msm, 0, 19, times, covariates)
expect_equal(0.0510873669808412, pp[1,3], tol=1e-04)
p <- pmatrix.msm(psor.msm, 5, covariates=covariates[[1]]) %*% pmatrix.msm(psor.msm, 5, covariates=covariates[[2]]) %*% pmatrix.msm(psor.msm, 5, covariates=covariates[[3]]) %*% pmatrix.msm(psor.msm, 4, covariates=covariates[[4]])
expect_equal(pp[1,3], p[1,3], tol=1e-04)
covariates <- list(list(ollwsdrt=0, hieffusn=0),list(ollwsdrt=0, hieffusn=1))
p <- pmatrix.piecewise.msm(psor.msm, 0, 7, times=5, covariates)
expect_equal(0.172773087945103, p[1,3], tol=1e-04)
expect_error(pmatrix.piecewise.msm("foo", 1,2, c(1, 2), c(1,2)), "expected .+ msm model")
expect_error(pmatrix.piecewise.msm("foo", 1, 2, c(1, 0.5, 2), list(0, 1, 0, 1)), "expected .+ msm model")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 0.5, 2), list(0, 1, 0, 1)), "times should be a vector of numbers in increasing order")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, "rubbish", list(0, 1, 0, 1)),"times should be a vector of numbers in increasing order")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 1.5, 2), "rubbish"))
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 1.5, 2), list("rubbish", "foo","bar","boing")), "covariates argument")
expect_error(pmatrix.piecewise.msm(psor.msm, 1, 2, c(1, 1.5, 2), list(0, 1, 0, 1)), "covariates argument")
})
test_that("pmatrix.piecewise.msm given just Q matrices",{
p1 <- pmatrix.piecewise.msm(
t1 = 0,
t2 = 6,
times = c(2, 4),
covariates = list(list(cov1 = 1), list(cov1 = 2), list(cov1 = 3)),
qlist = list(
Q1 = rbind(c(0.9, 0.1), c(0.1, 0.9)),
Q2 = rbind(c(0.8, 0.2), c(0.1, 0.9)),
Q3 = rbind(c(0.7, 0.3), c(0.1, 0.9))
)
)
expect_equal(round(p1), rbind(c(166, 238), c(99, 304)))
})
test_that("totlos.msm",{
tl <- totlos.msm(psor.msm, ci="normal", B=2)
sj <- sojourn.msm(psor.msm, ci="none")
expect_equal(tl[1,1], sj$estimates[[1]])
en <- envisits.msm(psor.msm, ci="normal", B=2)
expect_equal(en[1,2], 1)
})
test_that("logLik.msm",{
expect_equivalent(unclass(logLik.msm(psor.msm)), psor.msm$minus2loglik / -2)
expect_error(logLik.msm("foo"))
})
test_that("lrtest.msm",{
psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~ollwsdrt, constraint = list(ollwsdrt=c(1,1,2)), control=list(fnscale=1))
expect_equal(lrtest.msm(psor2.msm, psor.msm)[1,"-2 log LR"],
psor2.msm$minus2loglik - psor.msm$minus2loglik)
})
test_that("qmatrix subset function",{
Q <- qmatrix.msm(psor.msm)
expect_equivalent(Q[1,2]["SE"], Q$SE[1,2])
expect_error(Q[1,2,3,4,5], "unused arguments")
expect_error(Q[1], "Two dimensions must be supplied")
expect_equivalent(Q[], Q)
Q[1,]
Q[,2]
Q[c(1,2),]
Q[c(2,1),]
Q[c(1,2),c(1,2)]
})
test_that("efpt.msm",{
Q <- twoway4.q
expect_equal(efpt.msm(qmatrix=Q, tostate=3), c(Inf,Inf,0,Inf))
Q <- rbind(c(-0.25,0.25,0), c(0.166, -0.332, 0.166), c(0, 0.25, -0.25))
expect_equal(efpt.msm(qmatrix=Q, tostate=3)[c(1,2)], solve(-Q[1:2,1:2], c(1,1)),tol=1e-06)
Q <- twoway4.q; Q[2,4] <- Q[2,1] <- 0; diag(Q) <- 0; diag(Q) <- -rowSums(Q)
expect_equal(efpt.msm(qmatrix=Q, tostate=3), c(Inf, 6.02409638554217, 0, Inf), tol=1e-05)
expect_equal(efpt.msm(psor.msm, tostate=c(2)), c(10.4237243422138, 0, Inf, Inf), tol=1e-05)
expect_equal(efpt.msm(psor.msm, tostate=c(3)), c(16.5099104995137, 6.08618615729989, 0, Inf), tol=1e-05)
expect_equal(efpt.msm(psor.msm, tostate=c(2,3)), c(10.4237243422138, 0, 0, Inf), tol=1e-05)
})
test_that("ppass.msm",{
pp <- ppass.msm(psor.msm, tot=10)
pm <- pmatrix.msm(psor.msm, t=10)
expect_equal(pp[,4], pm[,4]) # state 4 is absorbing
pp <- ppass.msm(qmatrix=twoway4.q, tot=1000)
expect_equal(pp[1,2], 0.5)
expect_warning(ppass.msm(qmatrix=twoway4.q, tot=100, ci="normal"), "No fitted model supplied: not calculating confidence intervals")
})
test_that("score residuals",{
sres <- scoreresid.msm(psor.msm)
expect_equal(c(0.0608163009112361, 0.187998750251689, 0.0143302186951471), as.numeric(sres[1:3]), tol=1e-05)
psor2 <- na.omit(psor); psor2$months[psor2$ptnum==5] <- psor2$months[psor2$ptnum==5]*10
psor.0.q <- rbind(c(0,0.1,0,0),c(0,0,0.2,0),c(0,0,0,0.3),c(0,0,0,0))
psor.infl.msm <- msm(state ~ months, subject=ptnum, data=psor2, covariates = ~ollwsdrt+hieffusn, constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), qmatrix = psor.0.q)
sres <- scoreresid.msm(psor.infl.msm)
if (interactive()) sres <- scoreresid.msm(psor.infl.msm, plot=TRUE)
expect_equal(names(which.max(sres)), "5")
})
test_that("observed, expected etc",{
## two covariates
expect_equal(get.covhist(psor.msm)$example$time, c(6.4606, 17.078, 26.3217, 30.5763, 6.4052, 7.6893, 3.3593, 12.7775))
expect_equivalent(observed.msm(psor.msm)$obstab[,"State 4"], c(0, 6, 25, 39, 48, 50, 52, 53, 55, 56, 57))
expect_equal(as.numeric(expected.msm(psor.msm, covariates="population")$Expected[1:5,"State 3"]),
c(0, 11.7058415609289, 14.8078062584701, 9.50355791490463, 5.83372436779321), tol=1e-05)
expect_equal(as.numeric(expected.msm(psor.msm, covariates="mean")$Expected[1:5,"State 3"]),
c(0, 11.165395364321, 14.8732589026358, 9.67192721319246, 6.26754773418822), tol=1e-05)
})
test_that("observed, expected etc with one covariate, should run",{
expect_error({
psor1.msm <- msm(state ~ months, covariates=~ollwsdrt, subject=ptnum, data=psor, qmatrix = psor.q)
get.covhist(psor1.msm)
observed.msm(psor1.msm)
expected.msm(psor1.msm)
expected.msm(psor1.msm, covariates="mean")
}, NA)
})
test_that("observed, expected etc with PCI, should run",{
expect_error({
psor1.msm <- msm(state ~ months, covariates=~ollwsdrt+hieffusn, subject=ptnum, data=psor, qmatrix = psor.q, pci=c(5,10))
covhist <- get.covhist(psor1.msm)
observed.msm(psor1.msm)
expected.msm(psor1.msm)
expected.msm(psor1.msm, covariates="mean")
}, NA)
})
test_that("subset argument to observed",{
subs <- psor.msm$data$mf$"(subject)"[!duplicated(psor.msm$data$mf$"(subject)") & psor.msm$data$mf$ollwsdrt==0]
expect_equivalent(observed.msm(psor.msm, subset=subs)$obstab[,"State 4"], c(0, 5, 20, 30, 37, 38, 40, 41, 43, 44, 45))
})
test_that("error handling: formula",{
## formula
expect_error(msm(), "state ~ time formula not given")
expect_error(cav.msm <- msm(state, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "not found")
expect_error(cav.msm <- msm(~1, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "invalid data")
expect_error(cav.msm <- msm("foo", subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE), "not a formula")
})
test_that("error handling: qmatrix",{
wrong.q <- cbind(c(0,1,2), c(0,1,2))
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"Number of rows and columns of qmatrix should be equal")
wrong.q <- cbind(c(0,1), c(0,1))
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, fixedpars=TRUE),"State vector contains elements not in 1, 2")
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"Not all the states specified in \"deathexact\" are absorbing")
wrong.q <- "foo"
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"qmatrix should be a numeric matrix")
wrong.q <- 1
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, deathexact = TRUE, fixedpars=TRUE),"qmatrix should be a numeric matrix")
})
test_that("error handling: subject",{
expect_error(cav.msm <- msm(state~years, subject="foo", data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),
"variable lengths differ")
expect_error(cav.msm <- msm(state~years, subject=foo, data = cav,
qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),
"not found")
})
test_that("error handling: obstype",{
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype="foo", deathexact = TRUE, fixedpars=TRUE),"should be numeric")
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=rep(1,10), deathexact = TRUE, fixedpars=TRUE),"obstype of length") ##FIXME
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=rep(4, nrow(cav)), deathexact = TRUE, fixedpars=TRUE),"elements of obstype should be 1, 2, or 3")
obstype <- rep(1, nrow(cav))
obstype[c(1,8)] <- 5
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=obstype, deathexact = TRUE, fixedpars=TRUE), NA) # no error: obstype for first subject doesn't matter
obstype[2] <- 5
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=obstype, deathexact = TRUE, fixedpars=TRUE),"elements of obstype should be 1, 2, or 3") # error
})
test_that("error handling: covariates",{
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = "wibble"),"should be a formula or list of formulae")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sux),"not found")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=TRUE, misccovariates = "wobble"),"should be a formula")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sox),"not found")
})
test_that("error handling: covinits",{
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits="foo", fixedpars=TRUE),"covinits should be a list")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(sex="foo", age="bar"), fixedpars=TRUE),"should be numeric")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(sex=c(1,2,3), age="bar"), fixedpars=TRUE),"should be a list of numeric vectors")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(age=rep(0.1, 7)), fixedpars=TRUE),"covariate age in covinits unknown")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(age=rep(0.1, 7), foo=1, bar=2), fixedpars=TRUE),"covariates age, foo, bar in covinits unknown")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, covinits=list(sex=1), fixedpars=TRUE),"initial values of length 1, should be 7")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, misccovinits="foo", fixedpars=TRUE), "hcovinits should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, misccovariates = ~ sex, misccovinits=list(sex=1, age="bar"), fixedpars=TRUE), "misccovariates supplied but no ematrix")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, misccovariates = ~ sex, misccovinits=list(sex=1, age="bar"), fixedpars=TRUE),"misccovariates supplied but no ematrix")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, ematrix=ematrix, misccovariates = ~ sex, misccovinits=list(sex=1, age="bar"), fixedpars=TRUE),"covinits should be a list of numeric")
})
test_that("error handling: constraints",{
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, constraint="foo"),"constraint should be a list")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, constraint=list(foo="bar")),"constraint should be a list of numeric vectors")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, covariates = ~ sex, constraint=list(foo=1)),"Covariate .+ in constraint statement not in model.")
expect_warning(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, constraint="foo", fixedpars=TRUE),"constraint specified but no covariates")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, miscconstraint=list(foo="bar")),"constraint should be a list of numeric vectors")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, miscconstraint=list(foo=1)),"Covariate .+ in constraint statement not in model.")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, misccovariates = ~ sex, miscconstraint=list(sex=1)),"constraint of length 1, should be 4")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, qconstraint="foo", fixedpars=TRUE),"qconstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, qconstraint=list(c(1,1,2)), fixedpars=TRUE),"qconstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, qconstraint=c(1,1,2), fixedpars=TRUE),"baseline intensity constraint of length 3, should be 7")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, econstraint="foo", fixedpars=TRUE) ,"econstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, econstraint=list(c(1,1,2)), fixedpars=TRUE) ,"econstraint should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, econstraint=c(1,1,2), fixedpars=TRUE),"baseline misclassification constraint of length 3, should be 4")
})
test_that("error handling: initprobs",{
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, initprobs="poo", fixedpars=TRUE),"initprobs should be numeric")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, initprobs=c(1,2), fixedpars=TRUE),"initprobs vector of length 2, should be vector of length 4 or a matrix")
expect_error(cav.msm <- msm( state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.q, ematrix=ematrix, initprobs=c(2,1,1,1), fixedpars=TRUE), NA) # scaled to sum to 1.
})
test_that("error handling: check states",{
wrong.q <- cbind(c(0,1), c(0,1)) # extra states in data
expect_error(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, fixedpars=TRUE),"State vector contains elements not in 1, 2")
wrong.q <- rbind(c(0,1,2,3,1), c(0,1,3,4,1), c(0,1,2,3,2), c(0,1,2,3,4), c(0,0,0,0,0))
expect_warning(cav.msm <- msm(state~years, subject=PTNUM, data = cav, qmatrix = wrong.q, fixedpars=TRUE),"State vector doesn't contain observations of 5")
})
test_that("error handling: check times",{
cav.wrong <- cav
cav.wrong$years[3:5] <- 4:2
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),"not ordered by time")
cav.wrong <- cav
cav.wrong$PTNUM[4:5] <- 100003
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE),"Observations within subjects .+ are not adjacent in the data")
cav2 <- cav
cav2$years[10] <- cav2$years[9]
expect_warning(msm(state ~ years, subject=PTNUM, data = cav2, qmatrix = twoway4.q, fixedpars=TRUE),
"Different states observed at the same time on the same subject at observations 9 and 10")
## with missing data
cav2$years[6] <- NA ## report original rows before excluding missing data
expect_warning(msm(state ~ years, subject=PTNUM, data = cav2, qmatrix = twoway4.q, fixedpars=TRUE),
"Different states observed at the same time on the same subject at observations 9 and 10")
cav2 <- cav2[6:nrow(cav),]
expect_warning(msm(state ~ years, subject=PTNUM, data = cav2, qmatrix = twoway4.q, fixedpars=TRUE),
"Different states observed at the same time on the same subject at observations 4 and 5",
"Subject 100002 only has one complete observation")
})
test_that("error handling: check model",{
cav.wrong <- cav
cav.wrong$state[4] <- 1
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = oneway4.q, deathexact = TRUE, fixedpars=TRUE),"Data may be inconsistent with transition matrix for model without misclassification:\nindividual 100002 moves from state 2 to state 1 at observation 4")
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = oneway4.i, gen.inits=TRUE, fixedpars=TRUE),"individual 100046 moves from state 2 to state 1 at observation 225")
## row number reporting with missing data - should be the same
cav.wrong$PTNUM[2] <- NA
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = oneway4.q, deathexact = TRUE, fixedpars=TRUE),"100002 moves from state 2 to state 1 at observation 4")
cav.wrong <- cav
cav.wrong$state[4] <- 1
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, exacttimes=TRUE, fixedpars=TRUE),"individual 100003 moves from state 1 to state 3 at observation 10")
## row number reporting with missing data
cav.wrong$state[3] <- NA
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, exacttimes=TRUE, fixedpars=TRUE),"individual 100003 moves from state 1 to state 3 at observation 10") # should complain about obs 10
obstype <- rep(2, nrow(cav))
obstype[10] <- 1
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, obstype=obstype, fixedpars=TRUE),"individual 100006 moves from state 1 to state 3 at observation 29")
## absorbing-absorbing transitions
cav.wrong$state[6] <- 4
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.wrong, qmatrix = twoway4.q, obstype=obstype, fixedpars=TRUE),"Absorbing - absorbing transition at observation 7")
})
test_that("error handling: death",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = "foo", fixedpars=TRUE) ,"Exact death states indicator must be numeric")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = 5, fixedpars=TRUE) ,"Exact death states indicator contains states not in 1, 2, ... , 4")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = 1:5, fixedpars=TRUE) ,"Exact death states indicator contains states not in 1, 2, ... , 4")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, deathexact = 3, fixedpars=TRUE),"Not all the states specified in \"deathexact\" are absorbing" )
})
test_that("error handling: censor",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, censor="rubbish", fixedpars=TRUE),"censor must be numeric")
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, censor=1, fixedpars=TRUE),"some censoring indicators are the same as actual states")
expect_warning(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, censor.states="rubbish", fixedpars=TRUE) ,"censor.states supplied but censor not supplied")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.cens, qmatrix = twoway4.q, censor=99, censor.states="rubbish", fixedpars=TRUE) ,"censor.states should be all numeric")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav.cens, qmatrix = twoway4.q, censor=99, censor.states=list(c(1,2,3), "rubbish"), fixedpars=TRUE) ,"censor.states should be a vector")
})
test_that("error handling: obstype",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype="rubbish", fixedpars=TRUE) ,"obstype should be numeric")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=c(1,2,3), fixedpars=TRUE) ,"obstype of length 3, should be length 1 or 2846")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=4, fixedpars=TRUE),"elements of obstype should be 1, 2, or 3" )
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, obstype=rep(4,nrow(cav)), fixedpars=TRUE) ,"elements of obstype should be 1, 2, or 3")
})
test_that("error handling: fixedpars",{
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars="foo"),"Elements of fixedpars should be in 1, ..., 7")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=list(c(1,3,4))),"Elements of fixedpars should be in 1, ..., 7")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=1:8),"Elements of fixedpars should be in 1, ..., 7")
expect_error(cav.msm <- msm(state ~ years, subject=PTNUM, data = cav, qmatrix = twoway4.q, fixedpars=0.5),"Elements of fixedpars should be in 1, ..., 7")
})
test_that("error handling: plot",{
expect_error(plot.msm("foo"),"expected .+ msm model")
expect_error(plot.msm(psor.msm, from="foo"), "from must be numeric")
expect_error(plot.msm(psor.msm, to="foo"), "to must be numeric")
expect_error(plot.msm(psor.msm, from = 1:8, to=3),"from must be a vector of states in 1, ..., 4")
expect_error(plot.msm(psor.msm, to = 3),"to must be an absorbing state")
expect_error(plot.msm(psor.msm, range="foo"))
expect_error(plot.msm(psor.msm, range=1:6),"range must be a numeric vector of two elements")
})
test_that("recreate.olddata",{
expect_error(recreate.olddata(psor.msm), NA)
})
test_that("form.output",{
qo <- msm.form.qoutput(psor.msm)
expect_equal(qo$base.Estimate[2],
qmatrix.msm(psor.msm, covariates="mean")$estimates[1,2])
print(qmatrix.msm(psor.msm))
})
test_that("plots",{
skip_on_cran()
expect_error({
plot.msm(psor.msm)
plotprog.msm(state ~ months, subject=ptnum, data=psor)
plot.survfit.msm(psor.msm)
contour.msm(psor.msm)
persp.msm(psor.msm)
image.msm(psor.msm)
surface.msm(psor.msm, xrange=c(-2.4,-2.3))
}, NA)
})
test_that("miscellaneous",{
expect_error(msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~ollwsdrt+hieffusn,
constraint = list(hieffusn=c(1,1,1),ollwsdrt=c(1,1,2)), control=list(fnscale=1),
na.action=na.fail))
expect_error(model.matrix(psor.msm), NA)
expect_error(model.frame(psor.msm), NA)
})
test_that("single-column matrix in covariates",{
psor.q <- rbind(c(0,0.1,0,0),c(0,0,0.1,0),c(0,0,0,0.1),c(0,0,0,0))
psor$matcov <- matrix(psor$ollwsdrt, ncol=1)
psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~matcov, fixedpars=TRUE)
psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q,
covariates = ~ollwsdrt, fixedpars=TRUE)
expect_equal(psor.msm$minus2loglik, psor2.msm$minus2loglik)
})
Try the msm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.