-
Notifications
You must be signed in to change notification settings - Fork 17
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
6 changed files
with
188 additions
and
98 deletions.
There are no files selected for viewing
This file was deleted.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,34 @@ | ||
| 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)) | ||
|
|
||
| ## Merge states 2 and 3 | ||
| psor$state3 <- ifelse(psor$state==3, 2, | ||
| ifelse(psor$state==4, 3, psor$state)) | ||
| psor3.q <- psor.q[1:3,1:3] | ||
| psor3.msm <- msm(state3 ~ months, subject=ptnum, data=psor, qmatrix = psor3.q, | ||
| covariates = ~ollwsdrt+hieffusn, | ||
| constraint = list(hieffusn=c(1,1),ollwsdrt=c(1,1)), | ||
| control=list(fnscale=1)) | ||
|
|
||
| test_that("DRAIC",{ | ||
| d <- draic.msm(psor.msm, psor3.msm) | ||
| expect_true(is.numeric(d$draic)) | ||
| dl <- draic.msm(psor.msm, psor3.msm, likelihood.only=TRUE) | ||
| expect_equal(d$lik.restricted["complex","-LL"], dl[["complex"]]) | ||
| expect_equal(as.numeric(dl[["complex"]]), 415.032735368145, tolerance=1e-06) | ||
| expect_true(is.numeric(draic.msm(psor.msm, psor3.msm, tl=0.99)$ti["Lower"])) | ||
| }) | ||
|
|
||
| test_that("DRAIC with observed information", { | ||
| skip_on_cran() | ||
| expect_error(draic.msm(psor.msm, psor3.msm, information="observed"), NA) | ||
| }) | ||
|
|
||
| test_that("DRLCV",{ | ||
| psor.msm <- msm(state ~ months, subject=ptnum, data=psor[1:29,], qmatrix = psor.q) | ||
| psor3.msm <- msm(state3 ~ months, subject=ptnum, data=psor[1:29,], qmatrix = psor3.q) | ||
| expect_error(drlcv.msm(psor.msm, psor3.msm, verbose=FALSE), NA) | ||
| }) | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,14 @@ | ||
| test_that("pearson.msm help example",{ | ||
| 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.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))) | ||
| p <- pearson.msm(psor.msm, timegroups=2, intervalgroups=2, covgroups=2) | ||
| expect_equal(p$test[["p"]],0) | ||
| }) | ||
|
|
||
| #test_that("pearson with bootstrap",{ | ||
| # psor.msm <- msm(state ~ months, subject=ptnum, data=psor[1:29,], | ||
| # qmatrix = psor.q) | ||
| # p <- pearson.msm(psor.msm, timegroups=2, intervalgroups=2, covgroups=1, boot=TRUE, B=10) | ||
| #}) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,102 @@ | ||
| test_that("one phased state", { | ||
| psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, fixedpars=TRUE, | ||
| phase.states=c(1)) | ||
| expect_equal(psor.msm$minus2loglik, 1308.35590680014, tol=1e-06) | ||
| } | ||
| ) | ||
|
|
||
| test_that("working example",{ | ||
| skip_on_cran() | ||
| ## SIMULATE | ||
| mst1 <- 5 # Short-stay mean | ||
| mst2 <- 30 # Long-stay mean | ||
| p2 <- 0.9 # Long-stay probability | ||
| q23 <- 1/5 # Transition rate between phases | ||
|
|
||
| l1 <- (p2/mst1) | ||
| mu1 <- (1-p2)/mst1 | ||
| mu2 <- 1/(mst2-mst1) | ||
| Q <- rbind(c(0,l1,mu1*0.4,mu1*0.6), | ||
| c(0,0,mu2*0.4,mu2*0.6), | ||
| c(0,0,0,q23), | ||
| c(0,0,0,0)) | ||
| # Given the hidden state, the observed state is deterministic | ||
| E <- rbind(c(1,0,0,0), | ||
| c(1,0,0,0), | ||
| c(0,1,0,0), | ||
| c(0,0,1,0)) | ||
| nsubj <- 1000; nobspt <- 10 | ||
| set.seed(1) | ||
| sim.df <- data.frame(subject = rep(1:nsubj, each=nobspt), time = seq(0, 100, length=nobspt)) | ||
| sim2.df <- simmulti.msm(sim.df[,1:2], qmatrix=Q, ematrix=E) | ||
| statetable.msm(obs, subject, sim2.df) | ||
|
|
||
| ## FIT | ||
| Q3 <- rbind(c(0,0.5,0.5),c(0,0,0.5),c(0,0,0)) | ||
| s.msm <- msm(obs ~ time, subject=subject, data=sim2.df, phase.states=1, qmatrix=Q3, # default inits | ||
| phase.inits=list(list(trans=0.05, | ||
| exit=matrix(c(0.1,0.1,0.1,0.1),nrow=2))), | ||
| control=list(trace=1,REPORT=1,fnscale=50000,maxit=10000),method="BFGS") | ||
|
|
||
| # Parameter estimates should agree with the values used for simulation | ||
| s.msm | ||
| c(l1, mu1*0.4, mu1*0.6, mu2*0.4, mu2*0.6, q23) | ||
| means <- phasemeans.msm(s.msm) | ||
|
|
||
| expect_equal(means[[1,"Short stay mean"]], 5, tol=1) | ||
| expect_equal(means[[1,"Long stay mean"]], 30, tol=5) | ||
| expect_equal(means[[1,"Long stay probability"]], 0.9, tol=0.1) | ||
| }) | ||
|
|
||
|
|
||
| test_that("two phased states", { | ||
| ## note phase transition within state 2 appears hard to estimate. | ||
| psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, fixedpars=TRUE, | ||
| phase.states=c(1,2)) | ||
| expect_equal(psor.msm$minus2loglik, 1303.768173132893, tol=1e-06) | ||
| }) | ||
|
|
||
|
|
||
| test_that("supplying initial values",{ | ||
| psor.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1,2), | ||
| phase.inits = list( | ||
| list(trans=c(0.42032), exit=c(0.17404,0.04809)), | ||
| list(trans=c(0.42032), exit=c(0.17404,0.04809))), | ||
| fixedpars=TRUE) | ||
| expect_equal(psor.msm$minus2loglik, 1280.788860486683, tol=1e-06) | ||
| }) | ||
|
|
||
| test_that("errors in initial values",{ | ||
| expect_error(psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1,2), phase.inits = "foo"), "phase.inits should be a list") | ||
| expect_error(psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1,2), phase.inits = list(1,2,3)), "phase.inits of length 3, but there are 2 phased states") | ||
| expect_error(psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1,2), phase.inits = list(1,2)), "phase.inits.+ list of length 1, should be 2") | ||
| expect_error(psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1), phase.inits = list(list(trans=c(0.42032, 1), exit=c(0.17404,0.04809)))), "phase.inits.+trans of length 2, should be 1") | ||
| expect_error(psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1), phase.inits = list(list(trans=c(0.42032), exit=c(0.17404,0.04809,1,1)))), "phase.inits.+exit has 4 columns, should be 2") | ||
| expect_error(psor2.msm <- msm(state ~ months, subject=ptnum, data=psor, qmatrix = psor.q, phase.states=c(1), phase.inits = list(list(trans=c(0.42032), exit=matrix(c(0.17404,0.04809,1,1),ncol=2,byrow=TRUE)))), "phase.inits.+exit has 2 rows, but there are 1 exit states from this state") | ||
| }) | ||
|
|
||
|
|
||
| test_that("with censoring", { | ||
| psor2 <- psor | ||
| psor2$state[c(16,18)] <- 99 | ||
| psorc.msm <- msm(state ~ months, subject=ptnum, data=psor2, qmatrix = psor.q, fixedpars=TRUE, | ||
| censor = 99, censor.states=c(3,4)) | ||
| expect_equal(psorc.msm$minus2loglik, 1290.189703452163, tol=1e-06) | ||
| }) | ||
|
|
||
|
|
||
| test_that("HMM on top", { | ||
| ### Note 1.6 clarifies obstrue not supported for hmmCat: so test result in 1.5 wrong | ||
| miscnew.msm <- msm(state ~ years, subject = PTNUM, data = cav, | ||
| qmatrix = oneway4.q, death = 5, #obstrue=firstobs, | ||
| fixedpars=TRUE, | ||
| # could try to fit, but weakly identifiable. fit seems improved though (-2LL 3864 vs 3951) | ||
| # control=list(trace=1,REPORT=1,fnscale=4000, maxit=10000), | ||
| phase.states = 1, | ||
| hmodel=list( | ||
| hmmCat(prob=c(0.9, 0.1, 0, 0)), | ||
| hmmCat(prob=c(0.1, 0.8, 0.1, 0)), | ||
| hmmCat(prob=c(0, 0.1, 0.9, 0)), hmmIdent(4)) | ||
| ) | ||
| expect_equal(miscnew.msm$minus2loglik, 4357.70908245511, tol=1e-06) | ||
| }) |