cdata = data.frame(x2 = runif(nn <- 1000)) # ystar are true values
cdata = transform(cdata, ystar = rnorm(nn, m = 100 + 15 * x2, sd = exp(3)))
with(cdata, hist(ystar))
cdata = transform(cdata, L = runif(nn, 80, 90), # Lower censoring points
U = runif(nn, 130, 140)) # Upper censoring points
cdata = transform(cdata, y = pmax(L, ystar)) # Left censored
cdata = transform(cdata, y = pmin(U, y)) # Right censored
with(cdata, hist(y))
Extra = list(leftcensored = with(cdata, ystar < L),
rightcensored = with(cdata, ystar > U))
fit1 = vglm(y ~ x2, cennormal1, cdata, crit = "c", extra = Extra, trace = TRUE)
fit2 = vglm(y ~ x2, tobit(Lower = with(cdata, L), Upper = with(cdata, U)),
cdata, crit = "c", trace = TRUE)
coef(fit1, matrix = TRUE)
max(abs(coef(fit1, matrix = TRUE) - coef(fit2, matrix = TRUE))) # Should be 0
names(fit1@extra)
Run the code above in your browser using DataLab