library("truncSP") ########################################### ### plots ########################################### set.seed(19461203) ## Simulated data (for figures 1-4) n <- 200 x <- runif(n, 0, 10) eps <- rnorm(n, 0, 6) y <- 0 + 2 * x + eps xobs <- x[y > 0] yobs <- y[y > 0] xunobs <- x[y <= 0] yunobs <- y[y <= 0] ### Figure 1 ols <- lm(yobs ~ xobs) olsfull <- lm(y ~ x) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-15, 42), pch = 20) par(new = TRUE) plot(xunobs, yunobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-15, 42), pch = 4) abline(h = 0) abline(coef = ols$coef, lty = 2) abline(coef = olsfull$coef, lty = 3) text(6.5, -8, "Truncated observations", cex = 1.5) text(2, 13.2, expression(y == x^T * hat(beta)[obs]), cex = 1.5) text(4, 3.7, expression(y == x^T * hat(beta)[comp]), cex = 1.5) ## Figure 2 STLS par(mai = c(1.36, 1.093333, 1.093333, 1.56)) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0.3, 10), ylim = c(-15, 42), type = "n") polygon(c(0, 11, 11), c(0, 0, 44), border = NA, col = "grey80") polygon(c(-0.05, -0.05, 11, 11), c(0, 44, 44, 44), border = NA, col = "grey95") abline(h = 0) abline(h = 0, lty = 2, lwd = 2) abline(a = 0, b = 2) abline(a = 0, b = 4, lty = 2, lwd = 2) par(new = TRUE) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0.3, 10), ylim = c(-15, 42), pch = 20) par(new = TRUE) plot(xunobs, yunobs, xlab = "x", ylab = "y", xlim = c(0.3, 10), ylim = c(-15, 42), pch = 4) mtext(expression(y == x^T * beta), cex = 1.5, side = 4, las = 1, at = c(21.7)) mtext(expression(2 * x^T * beta), cex = 1.5, side = 4, las = 1, at = c(42)) text(3, 22, expression(italic((y/2)^2)), cex = 1.2) text(8, 5.5, expression(italic((y - x^T * beta)^2)), cex = 1.2) ## Figure 3(a) QME qmemod <- qme(yobs ~ xobs) par(mai = c(1.36, 1.093333, 1.093333, 1.56)) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-11, 32), type = "n") gr <- qmemod$cval$Value/2 gr2 <- 22 - qmemod$cval$Value gr3 <- 22 + qmemod$cval$Value polygon(c(-1, -1, gr, gr), c(0, (4 * gr), (4 * gr), 0), border = NA, col = "grey80") polygon(c(gr, gr, 11, 11), c(0, 0, 0, gr2), border = NA, col = "grey95") polygon(c(-1, -1, 11, 11, gr), c((4 * gr), 44, 44, gr3, (4 * gr)), border = NA, col = "grey95") polygon(c(gr, gr, 11, 11), c(0, (4 * gr), gr3, gr2), border = NA, col = "grey65") abline(h = 0) segments(gr, 2 * gr, 11, 22) segments(-1, 2 * gr, gr, 2 * gr) segments(-1, (4 * gr), gr, (4 * gr), lty = 2, lwd = 2) segments(-1, 0, gr, 0, lty = 2, lwd = 2) segments(gr, 0, 11, gr2, lty = 2, lwd = 2) segments(gr, (4 * gr), 11, (gr2 + 4 * gr), lty = 2, lwd = 2) segments(gr, 0, gr, (4 * gr), lty = 2, lwd = 2) par(new = TRUE) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-11, 32), pch = 20) par(new = TRUE) plot(xunobs, yunobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-11, 32), pch = 4) mtext(expression(y == x^T * beta), side = 4, cex = 1.2, las = 1, at = c(21.3)) mtext(expression(x^T * beta - c), side = 4, cex = 1.2, las = 1, at = c(15.1)) mtext(expression(x^T * beta + c), side = 4, cex = 1.2, las = 1, at = c(27)) text(1, 10, expression(italic((y - c)^2 - c^2)), cex = 1.2) text(3.1, 21, expression(italic((y - x^T * beta)^2 - c^2)), cex = 1.2) arrows(3.6, 19.8, 3.6, 10.4, length = 0.1) ## Figure 3(b) LT par(mai = c(1.36, 1.093333, 1.093333, 1.56)) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-11, 32), type = "n") gr4 <- 22 + 4 * gr polygon(c(-1, -1, gr, gr), c(0, (6 * gr), (6 * gr), 0), border = NA, col = "grey80") polygon(c(gr, gr, 11, 11), c(0, 0, 0, gr2), border = NA, col = "grey95") polygon(c(-1, -1, 11, 11, gr), c((6 * gr), 44, 44, gr4, (6 * gr)), border = NA, col = "grey95") polygon(c(gr, gr, 11, 11), c(0, (6 * gr), gr4, gr2), border = NA, col = "grey65") abline(h = 0) segments(gr, 2 * gr, 11, 22) segments(-1, 2 * gr, gr, 2 * gr) segments(-1, (6 * gr), gr, (6 * gr), lty = 2, lwd = 2) segments(-1, 0, gr, 0, lty = 2, lwd = 2) segments(gr, 0, 11, gr2, lty = 2, lwd = 2) segments(gr, 0, gr, (6 * gr), lty = 2, lwd = 2) segments(gr, (6 * gr), 11, (gr2 + 6 * gr), lty = 2, lwd = 2) par(new = TRUE) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-11, 32), pch = 20) par(new = TRUE) plot(xunobs, yunobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-11, 32), pch = 4) mtext(expression(y == x^T * beta), cex = 1.2, at = 21.5, side = 4, las = 1) mtext(expression(x^T * beta - c[L]), cex = 1.2, at = 15.5, side = 4, las = 1) mtext(expression(x^T * beta + c[U]), cex = 1.2, at = 32.5, side = 4, las = 1) text(1, 12.5, expression(italic((y - c[L])^2/2)), cex = 1.2) text(6.5, 30, expression(italic((y - x^T * beta)^2/2)), cex = 1.2) text(3.160156, 26, expression(italic(c[U]^2/2)), cex = 1.2) text(7.5, 3, expression(italic(c[L]^2/2)), cex = 1.2) arrows(6.6, 28.7, 6.6, 20.5, length = 0.1) ## Figure 4 stlsmod <- stls(yobs ~ xobs) ltmod <- lt(yobs ~ xobs) plot(xobs, yobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-15, 42), pch = 20) par(new = TRUE) plot(xunobs, yunobs, xlab = "x", ylab = "y", xlim = c(0, 10), ylim = c(-15, 42), pch = 4) abline(h = 0) abline(coef = ols$coef) abline(coef = stlsmod$coef, lty = 2) abline(coef = qmemod$coef, lty = 3) abline(coef = ltmod$coef, lty = 5) text(0.5, 37.5, "OLS", cex = 0.8) text(0.5, 35, "STLS", cex = 0.8) text(0.5, 32.5, "QME", cex = 0.8) text(0.5, 30, "LT", cex = 0.8) segments(1, 37.5, 2, 37.5, lty = 1) segments(1, 35, 2, 35, lty = 2) segments(1, 32.5, 2, 32.5, lty = 3) segments(1, 30, 2, 30, lty = 5) ## Figure 5 library(lattice) data("PM10", package = "truncSP") PM10long <- reshape(PM10, direction = "long", varying = list(2:ncol(PM10)), times = colnames(PM10)[-1], v.names = "Value", timevar = "Variable") PM10long$Variable <- factor(PM10long$Variable, colnames(PM10)[-1]) lattice::xyplot(PM10 ~ Value | Variable, data = PM10long, xlab = "", scales = list(x = list(relation = "free"), y = list(alternating = 3)), par.settings = standard.theme("pdf", color = FALSE), as.table = TRUE, layout = c(3, 3)) ########################################### ### illustrative example ########################################### data("PM10", package = "truncSP") data("PM10trunc", package = "truncSP") ## Fitted models (all variables, Table 1) olsfull.all <- lm(PM10 ~ cars + temp + wind.speed + temp.diff + wind.dir + hour + day, data = PM10) olstrunc.all <- lm(PM10 ~ cars + temp + wind.speed + temp.diff + wind.dir + hour + day, data = PM10trunc) stls.all <- stls(PM10 ~ cars + temp + wind.speed + temp.diff + wind.dir + hour + day, data = PM10trunc, point = 2, control = list(maxit = 3000), covar = TRUE) qme.all <- qme(PM10 ~ cars + temp + wind.speed + temp.diff + wind.dir + hour + day, data = PM10trunc, point = 2, covar = TRUE, control = list(maxit = 4000)) lt.all <- lt(PM10 ~ cars + temp + wind.speed + temp.diff + wind.dir + hour + day, data = PM10trunc, point = 2, covar = TRUE, control = list(maxit = 2500)) summary(olsfull.all) summary(olstrunc.all) summary(stls.all) summary(qme.all) summary(lt.all) ## Fitted models (only sign. variables, Table 2) olsfull.sign <- lm(PM10 ~ cars + wind.speed, data = PM10) olstrunc.sign <- lm(PM10 ~ cars + wind.speed, data = PM10trunc) stls.sign <- stls(PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, covar = TRUE) qme.sign <- qme(PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, covar = TRUE) lt.sign <- lt(PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, covar = TRUE) summary(olsfull.sign) summary(olstrunc.sign) summary(stls.sign) summary(qme.sign) summary(lt.sign) ## Fitted models with different threshold values (Table 3) qme.0.5 <- qme(PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, const = 0.5, covar = TRUE) qme.1 <- qme.sign qme.2 <- qme(PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, const = 2, covar = TRUE) summary(qme.0.5) summary(qme.1) summary(qme.2) ## Output examples library("truncSP") data("PM10trunc", package = "truncSP") qmeobj <- qme(formula = PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, covar = TRUE) summary(qmeobj) qmeobj <- qme(formula = PM10 ~ cars + wind.speed, data = PM10trunc, point = 2, covar = FALSE) summary(qmeobj) vcov(qmeobj)