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_equivalent(psor.msm$Qmatrices$baseline, 
                      unclass(qmatrix.msm(psor.msm, covariates="mean", ci="none")))
    expect_equivalent(psor.nocen.msm$Qmatrices$baseline, 
                      unclass(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_equivalent(psor.nocen.msm$Qmatrices$baseline, 
                      unclass(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_equivalent(unclass(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_equivalent(unclass(efpt.msm(qmatrix=Q, tostate=3)), 
                      c(Inf, 6.02409638554217, 0, Inf), tol=1e-05)
    expect_equivalent(unclass(efpt.msm(psor.msm, tostate=c(2))),
                      c(10.4237243422138, 0, Inf, Inf), tol=1e-05)
    expect_equivalent(unclass(efpt.msm(psor.msm, tostate=c(3))), 
                      c(16.5099104995137, 6.08618615729989, 0, Inf), tol=1e-05)
    expect_equivalent(unclass(efpt.msm(psor.msm, tostate=c(2,3))),
                      c(10.4237243422138, 0, 0, Inf), tol=1e-05)
})

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)
})