Fixes for weights and bootstrapping

NAMESPACE

@@ -28,6 +28,7 @@ S3method(tidy,msm.prevalence)
 export(MatrixExp)
 export(absorbing.msm)
 export(boot.msm)
+export(coef.msm)
 export(crudeinits.msm)
 export(d2phase)
 export(deltamethod)
@@ -61,6 +62,7 @@ export(hmmTNorm)
 export(hmmUnif)
 export(hmmWeibull)
 export(hmodel2list)
+export(logLik.msm)
 export(lrtest.msm)
 export(msm)
 export(msm.form.eoutput)

R/msm.R

@@ -939,7 +939,12 @@ msm <- function(formula, subject=NULL, data=list(), qmatrix, gen.inits=FALSE,
     mf <- eval(temp, parent.frame())
 
     ## remember user-specified names for later (e.g. bootstrap/cross validation)
-    usernames <- c(state=all.vars(formula[[2]]), time=all.vars(formula[[3]]), subject=as.character(temp$subject), obstype=as.character(substitute(obstype)), obstrue=as.character(temp$obstrue))
+    usernames <- c(state=all.vars(formula[[2]]),
+                   time=all.vars(formula[[3]]),
+                   subject=as.character(temp$subject),
+                   subject.weights=as.character(temp$subject.weights),
+                   obstype=as.character(substitute(obstype)),
+                   obstrue=as.character(temp$obstrue))
     attr(mf, "usernames") <- usernames
 
     ## handle matrices in state outcome constructed in formula with cbind()

R/outputs.R

@@ -1405,6 +1405,7 @@ pnext.msm <- function(x, covariates="mean", ci=c("normal","bootstrap","delta","n
 #' @author C. H. Jackson \email{chris.jackson@@mrc-bsu.cam.ac.uk}
 #' @seealso \code{\link{msm}}
 #' @keywords models
+#' @export coef.msm
 #' @export
 coef.msm <- function(object, ...)
 {
@@ -1435,6 +1436,7 @@ coef.msm <- function(object, ...)
 #' @author C. H. Jackson \email{chris.jackson@@mrc-bsu.cam.ac.uk}
 #' @seealso \code{\link{msm}},\code{\link{lrtest.msm}}.
 #' @keywords models
+#' @export logLik.msm
 #' @export
 logLik.msm <- function(object, by.subject=FALSE, ...)
 {

R/pearson.R

@@ -340,7 +340,7 @@ pearson.msm <- function(x, transitions=NULL, timegroups=3, intervalgroups=3, cov
         qmatindex <- match(pastedc,pastedu)
         qmat <- array(0,dim=c(nst,nst,nouniq))
         for (i in 1:nouniq)
-            qmat[,,i] <- qmatrix.msm(x, covariates=as.list(uniq[i,]), ci="none")
+            qmat[,,i] <- qmatrix.msm(x, covariates=as.list(uniq[i,,drop=FALSE]), ci="none")
     }else{
         qmatindex <- rep(1,n)
         nouniq <- 1
@@ -744,6 +744,7 @@ pearson.boot.msm <- function(x, imp.times=NULL, transitions=NULL, timegroups=4,
   x$call$hessian <- x$call$death <- FALSE
   x$call$obstype <- NULL
   x$call$subject <- substitute(subject)
+  x$call$subject.weights <- substitute(subject.weights)
   i <- 1
   while (i <= B) {
     if (!is.null(imp.times))

R/simul.R

@@ -550,6 +550,7 @@ simfitted.msm <- function(x, drop.absorb=TRUE, drop.pci.imp=TRUE){
     } else misccov.effs <- NULL
     names(sim.df) <- replace(names(sim.df), match(c("(time)","(subject)"), names(sim.df)),
                              c("time","subject"))
+    if (!is.null(sim.df$"(subject.weights)")) names(sim.df)[names(sim.df)=="(subject.weights)"] = "subject.weights"
     sim.df$state <- NULL # replace observed with simulated state
     if (any(union(names(cov.effs), names(misccov.effs)) %in% c("state","time","subject","cens")))
         stop("Not supported with covariates named \"state\", \"time\", \"subject\" or \"cens\"") # TODO?

tests/testthat/test_weights.R

@@ -1,4 +1,3 @@
-
 test_that("subject weights",{
   cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
                   qmatrix = twoway4.q, deathexact = TRUE, fixedpars=TRUE)
@@ -41,22 +40,27 @@ test_that("subject weights",{
 
 test_that("subject weights, model fitting",{
   skip_on_cran()
-
   cav.msm <- msm( state ~ years, subject=PTNUM, data = cav,
                      qmatrix = twoway4.q, deathexact = TRUE, fixedpars=FALSE)
-  
+
   cav$swt <- 1
   cavwt1.msm <- msm( state ~ years, subject=PTNUM, data = cav,
-                     qmatrix = twoway4.q, deathexact = TRUE, fixedpars=FALSE,
+                     # work around boot not finding this in full test()
+                     qmatrix = rbind(c(-0.5,0.25,0,0.25),c(0.166, -0.498, 0.166, 0.166), 
+                                     c(0, 0.25, -0.5, 0.25), c(0, 0, 0, 0)),  
+
+                     deathexact = TRUE, fixedpars=FALSE,
                      subject.weights = swt)
   expect_equal(cav.msm$minus2loglik, cavwt1.msm$minus2loglik)
   pmatrix.msm(cavwt1.msm, ci="normal")
-
+  set.seed(1)
+  pmatrix.msm(cavwt1.msm, ci="boot", B=3)
+  pearson.msm(cavwt1.msm, boot = TRUE, B=2)
+
   cav$swt[cav$PTNUM <= 100020] <- 1.2
   cavwt2.msm <- msm( state ~ years, subject=PTNUM, data = cav,
                      qmatrix = twoway4.q, deathexact = TRUE, fixedpars=FALSE,
                      subject.weights = swt)
-  pmatrix.msm(cavwt2.msm, ci="normal")
-
+
   expect_true(cav.msm$minus2loglik != cavwt2.msm$minus2loglik)
 })