Open yushuf opened 6 years ago
The problem is a negative variance. This seems to happen only if who-infects-whom is not observed. How often does this problem occur (e.g., per thousand simulations)?
The problem appears when who-infected-whom is not observed with contact tracing true for external intercept and external log(shape) parameters interval estimates. The probabilities of happening the event are
P(xintercept_unobs_tr_w) = 0.497 and P(lxshape_unobs_tr_w) = 0.454
This is based on one thousand simulation for Wald CI. I never tried LR, in these cases LR stops running the program.
When there are no external infections, there should be no external intercept or shape parameters. Can you explain the second case where the problem appears again?
Sorry I made a mistake in description, That never happened when we ignored the external source of infection. I have updated my previous comment. Thanks for figure that out.
This appears to be numerical instability caused by the rapidly decreasing hazard of external infection, not an error in the code. Please re-run the simulations with exponential external contact intervals.
#Recently I was trying transreg with dist = "weibull", xdist = "weibull" options. I get the following errors while estimating the confidence intervals. First one is Wald and second one is LR.
**My message, this is not the only problem with LR which might happen with Wald. I tried two simulation and both of the times I get that error with Wald though in some different command line.
This is for your information. The detail code is as follows along with data**
data.zip
library(transtat) pdata <- read.csv("pdata.csv") xdata <- read.csv("xdata.csv") coef <- read.csv("coef.csv", header = TRUE)
pairwise data to be combined with external data
pdat <- pdata pdat$interval <- pdat$exit - pdat$entry pdat$external <- 0 pdat$trace <- 1 pdat$start <- 0 pdat$stop <- pdat$interval
external data
xdat <- xdata names(xdat)[names(xdat) == "x"] <- "susx"
xdat$external <- 1 xdat$infectious <- NA xdat$infx <- 0 xdat$interval <- xdat$exit # entry is at time 0 xdat$start <- xdat$entry xdat$stop <- xdat$exit
dat <- rbind(pdat, xdat)
AFT model all data
Who-infects-whom observed
obs <- transreg( Surv(interval, infector) ~ infx + susx + ext(external), sus = "susceptible", data = dat, dist = "weibull", xdist = "weibull" )
Who-infects-whom not observed
unobs <- transreg( Surv(interval, infset) ~ infx + susx + ext(external), sus = "susceptible", data = dat, dist = "weibull", xdist = "weibull" )
AFT model with contact tracing and left-truncation (delayed entry)
Who-infects-whom observed
obs_tr <- transreg( Surv(start, stop, infector) ~ infx + susx + ext(external), sus = "susceptible", data = dat, subset = (trace == 1), dist = "weibull", xdist = "weibull" )
Who-infects-whom not observed
unobs_tr <- transreg( Surv(start, stop, infset) ~ infx + susx + ext(external), sus = "susceptible", data = dat, subset = (trace == 1), dist = "weibull", xdist = "weibull" )
AFT model with no external data
Who-infects-whom observed
obs_nx <- transreg( Surv(interval, infector) ~ infx + susx, sus = "susceptible", data = dat, subset = (external == 0), dist = "weibull", xdist = "weibull" )
Who-infects-whom not observed
unobs_nx <- transreg( Surv(interval, infset) ~ infx + susx, sus = "susceptible", data = dat, subset = (external == 0), dist = "weibull", xdist = "weibull" )
AFT model with contact tracing but no left truncation
Who-infects-whom observed
obs_trx <- transreg( Surv(interval, infector) ~ infx + susx + ext(external), sus = "susceptible", dat = dat, subset = (trace == 1 | infector == 1), dist = "weibull", xdist = "weibull" )
Who-infects-whom not observed
unobs_trx <- transreg( Surv(interval, infset) ~ infx + susx + ext(external), sus = "susceptible", data = dat, subset = (trace == 1 | infector == 1), dist = "weibull", xdist = "weibull" )
ci_obs_w <-confint(obs) ci_obs_lr <- confint(obs, type = "lr") ci_unobs_w <-confint(unobs) ci_unobs_lr <- confint(unobs, type = "lr") ci_obs_tr_w <- confint(obs_tr) ci_obs_tr_lr <- confint(obs_tr, type = 'lr') ci_unobs_tr_w <- confint(unobs_tr) ci_unobs_tr_lr <- confint(unobs_tr, type = "lr")# ci_obs_nx_w <- confint(obs_nx) ci_obs_nx_lr <- confint(obs_nx, type = "lr") ci_unobs_nx_w <- confint(unobs_nx) ci_unobs_nx_lr <- confint(unobs_nx, type = "lr") ci_obs_trx_w <- confint(obs_trx) ci_obs_trx_lr <- confint(obs_trx, type = "lr") ci_unobs_trx_w <- confint(unobs_trx) ci_unobs_trx_lr <- confint(unobs_trx, type = "lr")