# Here, fit1 is a standard Tobit model and fit2 is a nonstandard Tobit model
Lower = 1; Upper = 4; set.seed(1)  # For the nonstandard Tobit model
tdata = data.frame(x2 = seq(-1, 1, len = (nn <- 100)))
meanfun1 = function(x) 0 + 2*x
meanfun2 = function(x) 2 + 2*x
tdata = transform(tdata,
  y1 = rtobit(nn, mean = meanfun1(x2)),  # Standard Tobit model 
  y2 = rtobit(nn, mean = meanfun2(x2), Lower = Lower, Upper = Upper))
with(tdata, table(y1 == 0)) # How many censored values?
with(tdata, table(y2 == Lower | y2 == Upper)) # How many censored values?
with(tdata, table(attr(y2, "cenL")))
with(tdata, table(attr(y2, "cenU")))
fit1 = vglm(y1 ~ x2, tobit, tdata, trace = TRUE,
            crit = "coeff")  # crit = "coeff" is recommended
coef(fit1, matrix = TRUE)
summary(fit1)
fit2 = vglm(y2 ~ x2, tobit(Lower = Lower, Upper = Upper, type.f = "cens"),
            tdata, crit = "coeff", trace = TRUE) # ditto
table(fit2@extra$censoredL)
table(fit2@extra$censoredU)
coef(fit2, matrix = TRUE)
# Plot the results
par(mfrow = c(2, 1))
plot(y1 ~ x2, tdata, las = 1, main = "Standard Tobit model",
     col = as.numeric(attr(y1, "cenL")) + 3,
     pch = as.numeric(attr(y1, "cenL")) + 1)
legend(x = "topleft", leg = c("censored", "uncensored"),
       pch = c(2, 1), col = c("blue", "green"))
legend(-1.0, 2.5, c("Truth", "Estimate", "Naive"),
       col = c("purple", "orange", "black"), lwd = 2, lty = c(1, 2, 2))
lines(meanfun1(x2) ~ x2, tdata, col = "purple", lwd = 2)
lines(fitted(fit1) ~ x2, tdata, col = "orange", lwd = 2, lty = 2)
lines(fitted(lm(y1 ~ x2, tdata)) ~ x2, tdata, col = "black",
      lty = 2, lwd = 2) # This is simplest but wrong!
plot(y2 ~ x2, tdata, las = 1, main = "Tobit model",
     col = as.numeric(attr(y2, "cenL")) + 3 +
           as.numeric(attr(y2, "cenU")),
     pch = as.numeric(attr(y2, "cenL")) + 1 +
           as.numeric(attr(y2, "cenU")))
legend(x = "topleft", leg = c("censored", "uncensored"),
       pch = c(2, 1), col = c("blue", "green"))
legend(-1.0, 3.5, c("Truth", "Estimate", "Naive"),
       col = c("purple", "orange", "black"), lwd = 2, lty = c(1, 2, 2))
lines(meanfun2(x2) ~ x2, tdata, col = "purple", lwd = 2)
lines(fitted(fit2) ~ x2, tdata, col = "orange", lwd = 2, lty = 2)
lines(fitted(lm(y2 ~ x2, tdata)) ~ x2, tdata, col = "black",
      lty = 2, lwd = 2) # This is simplest but wrong!Run the code above in your browser using DataLab