1 2sm.discontinuity <- function(x, y, h, hd, ...) { 3 4# A test for the presence of one or more discontinuities 5 6 x.name <- deparse(substitute(x)) 7 if (isMatrix(x)) x.names <- dimnames(x)[[2]] 8 y.name <- deparse(substitute(y)) 9 10 opt <- sm.options(list(...)) 11 data <- sm.check.data(x = x, y = y, ...) 12 x <- data$x 13 y <- data$y 14 n <- data$nobs 15 ndim <- data$ndim 16 opt <- data$options 17 if (ndim > 2) x <- x[, 1:2] 18 19 replace.na(opt, display, "lines") 20 replace.na(opt, se, TRUE) 21 replace.na(opt, band, TRUE) 22 replace.na(opt, test, TRUE) 23 replace.na(opt, col, "black") 24 replace.na(opt, df, 5) 25 if (ndim == 1) { 26 replace.na(opt, ngrid, 100) 27 replace.na(opt, xlab, x.name) 28 replace.na(opt, ylab, y.name) 29 replace.na(opt, xlim, range(x)) 30 replace.na(opt, ylim, range(y)) 31 if (length(opt$lty) == 1) 32 opt$lty <- c(opt$lty, opt$lty + 1) 33 } 34 else { 35 replace.na(opt, ngrid, 21) 36 dimn <- x.names 37 name.comp<-if(!is.null(dimn) & !all(dimn == "")) dimn 38 else {if (!is.null(attributes(x)$names)) attributes(x)$names 39 else outer(x.name, c("[1]", "[2]"), paste, sep = "")} 40 replace.na(opt, xlab, name.comp[1]) 41 replace.na(opt, ylab, name.comp[2]) 42 replace.na(opt, xlim, range(x[, 1])) 43 replace.na(opt, ylim, range(x[, 2])) 44 } 45 46 if(missing(h)) 47 h <- h.select(x, y, ...) 48 doublesmooth <- TRUE 49 if (missing(hd)) { 50 if (ndim == 1) { 51 hd <- h * sqrt(0.25) 52 h <- h * sqrt(0.75) 53 } 54 else { 55 hd <- h * sqrt(0.5) 56 h <- h * sqrt(0.5) 57 } 58 } 59 else if (all(hd == rep(0, ndim))) 60 doublesmooth <- FALSE 61 62 if (ndim == 1) 63 result <- sm.discontinuity.1d(x, y, h, hd, doublesmooth, opt) 64 else 65 result <- sm.discontinuity.2d(x, y, h, hd, doublesmooth, opt) 66 67 result$h <- h 68 if (doublesmooth) result$hd <- hd 69 invisible(result) 70 } 71 72sm.discontinuity.1d <- function(x, y, h, hd, doublesmooth, opt) { 73 74 y <- y[order(x)] 75 x <- sort(x) 76 n <- length(x) 77 78# Define z to be the mid-points of distinct x values. 79# Restrict z so that there are at least 5 observations, 80# over more than one design point, on each size of 81# every value. 82 83 z <- x[c(1, diff(x)) > 0] 84 nz <- length(z) 85 z <- (z[1:(nz-1)] + z[2:nz]) / 2 86 nz <- length(z) 87 flag <- rep(T, nz) 88 for (i in 1:nz) { 89 left <- x[x < z[i]] 90 right <- x[x > z[i]] 91 flag[i] <- (length(left) > 5 92 & length(right) > 5 93 & length(diff(left)[diff(left) > 0]) > 1 94 & length(diff(right)[diff(right) > 0]) > 1) 95 } 96 z <- z[flag] 97 nz <- length(z) 98 99 ghat.left <- vector("numeric", length = nz) 100 ghat.right <- vector("numeric", length = nz) 101 102 wd <- matrix(rep(z, rep(n, nz)), ncol = n, byrow = T) 103 wd <- wd - matrix(rep(x, nz), ncol = n, byrow = T) 104 w <- exp(-.5 * (wd / h)^2) 105 106 wl <- w * (sign(wd) + 1) / 2 107 s0 <- wl %*% rep(1, n) 108 s1 <- (wl * wd) %*% rep(1, n) 109 s2 <- (wl * wd^2 ) %*% rep(1, n) 110 wl <- wl * (matrix(rep(s2, n), ncol = n) - wd * matrix(rep(s1, n), ncol = n)) 111 wl <- wl / (matrix(rep(s2, n), ncol = n) * matrix(rep(s0, n), ncol = n) 112 - matrix(rep(s1, n), ncol = n)^2) 113 ghat.left <- wl %*% y 114 115 wr <- w * (1 - sign(wd)) / 2 116 s0 <- wr %*% rep(1, n) 117 s1 <- (wr * wd) %*% rep(1, n) 118 s2 <- (wr * wd^2 ) %*% rep(1, n) 119 wr <- wr * (matrix(rep(s2, n), ncol = n) - wd * matrix(rep(s1, n), ncol = n)) 120 wr <- wr / (matrix(rep(s2, n), ncol = n) * matrix(rep(s0, n), ncol = n) 121 - matrix(rep(s1, n), ncol = n)^2) 122 ghat.right <- wr %*% y 123 124 A <- sm.sigweight(x, rep(1, length(x))) / (n - 2) 125 w <- wl - wr 126 if (doublesmooth) { 127 ws <- sm.weight(z, z, hd) 128 w <- ws %*% w 129 ghat.left <- as.vector(ws %*% as.vector(ghat.left)) 130 ghat.right <- as.vector(ws %*% as.vector(ghat.right)) 131 } 132 133 shat <- sqrt(as.vector(t(as.matrix(y)) %*% A %*% as.matrix(y))) 134 s.e. <- as.vector(shat * sqrt((w^2) %*% rep(1, n))) 135 ts <- sum(((ghat.left - ghat.right) / s.e.) ^2 ) 136 A <- t(w) %*% diag((shat / s.e.)^2) %*% w - A * ts 137 p <- p.quad.moment(A, diag(n), 0, 0) 138 139 if (opt$display != "none") { 140 if (!opt$add) 141 plot(x, y, xlab = opt$xlab, ylab = opt$ylab, xlim = opt$xlim, ylim = opt$ylim, type = "n") 142 av <- (ghat.left + ghat.right) / 2 143 if (opt$band & opt$se) 144 polygon(c(z, rev(z)), c(av + s.e., rev(av - s.e.)), col = opt$col.band, border = FALSE) 145 lines(z, ghat.left, lty = opt$lty[1]) 146 lines(z, ghat.right, lty = opt$lty[2]) 147 points(x, y, col = opt$col.points, pch = opt$pch) 148 } 149 150 if (opt$verbose > 0) 151 cat("Test of continuity: significance = ", round(p, 3), "\n") 152 st.diff <- (ghat.left - ghat.right)/ s.e. 153 diffmat <- cbind(z, round(st.diff, 2))[abs(st.diff) > 2.5,] 154 # The following line forces a matrix when there is only one row in diffmat. 155 if (!is.matrix(diffmat)) diffmat <- matrix(diffmat, ncol = 2) 156 if ((opt$verbose > 0) & (nrow(diffmat) > 0)) { 157 cat("location st.diff\n") 158 for (i in 1:nrow(diffmat)) cat(diffmat[i, ], "\n") 159 } 160 161 invisible(list(p = p, sigma = shat, eval.points = z, st.diff = st.diff, diffmat = diffmat)) 162 163 } 164 165sm.discontinuity.2d <- function(x, y, h, hd, doublesmooth, opt, 166 nangles = 4, trim = 1/6, hull = FALSE) { 167 168 # Discontinuity detection with two covariates 169 170 n <- nrow(x) 171 del1 <- diff(range(x[,1])) * trim 172 del2 <- diff(range(x[,2])) * trim 173 x1grid <- seq(min(x[,1]) + del1, max(x[,1]) - del1, length = opt$ngrid) 174 x2grid <- seq(min(x[,2]) + del2, max(x[,2]) - del2, length = opt$ngrid) 175 ev.points <- cbind(x1grid, x2grid) 176 replace.na(opt, eval.points, ev.points) 177 eval.points <- opt$eval.points 178 ngrid <- nrow(eval.points) 179 weights <- rep(1, n) 180 181 wd1 <- matrix(rep(eval.points[, 1], n), ncol = n) 182 wd1 <- wd1 - matrix(rep(x[, 1], ngrid), ncol = n, byrow = TRUE) 183 wd2 <- matrix(rep(eval.points[, 2], n), ncol = n) 184 wd2 <- wd2 - matrix(rep(x[, 2], ngrid), ncol = n, byrow = TRUE) 185 w1 <- exp(-0.5 * (wd1 / h[1])^2) 186 w1 <- w1 * matrix(rep(weights, ngrid), ncol = n, byrow = TRUE) 187 w2 <- exp(-0.5 * (wd2 / h[2])^2) 188 wy <- matrix(rep(y, ngrid), ncol = n, byrow=TRUE) 189 190 a11 <- w1 %*% t(w2) 191 a12 <- (w1 * wd1) %*% t(w2) 192 a13 <- w1 %*% t(w2 * wd2) 193 a22 <- (w1 * wd1^2) %*% t(w2) 194 a23 <- (w1 * wd1) %*% t(w2 * wd2) 195 a33 <- w1 %*% t(w2 * wd2^2) 196 197 c1 <- w1 %*% t(w2 * wy) 198 c2 <- (w1 * wd1) %*% t(w2 * wy) 199 c3 <- w1 %*% t(w2 * wy * wd2) 200 201 beta1 <- sm.regression.invert(a22,a12,a23,a11,a13,a33,c2,c1,c3) 202 beta2 <- sm.regression.invert(a33,a23,a13,a22,a12,a11,c3,c2,c1) 203 204 wmask <- matrix(1, nrow = ngrid, ncol = ngrid) 205 206 if (hull) { 207 hull.points <- x[order(x[,1], x[,2]),] 208 dh <- diff(hull.points) 209 hull.points <- hull.points[c(TRUE, !((dh[,1] == 0) & (dh[,2] == 0))),] 210 hull.points <- hull.points[chull(hull.points), ] 211 nh <- nrow(hull.points) 212 gstep <- matrix(rep(eval.points[2, ] - eval.points[1,], nh), 213 ncol = 2, byrow = TRUE) 214 hp.start <- matrix(rep(eval.points[1, ], nh), ncol = 2, byrow = TRUE) 215 hull.points <- hp.start + gstep * round((hull.points - hp.start) / gstep) 216 hull.points <- hull.points[chull(hull.points), ] 217 grid.points <- cbind(rep(eval.points[, 1], ngrid), 218 rep(eval.points[, 2], rep(ngrid, ngrid))) 219 D <- diff(rbind(hull.points, hull.points[1, ])) 220 temp <- D[, 1] 221 D[,1] <- D[, 2] 222 D[,2] <- -temp 223 C <- as.vector((hull.points * D) %*% rep(1, 2)) 224 C <- matrix(rep(C, ngrid^2), nrow = ngrid^2, byrow = TRUE) 225 D <- t(D) 226 wmask <- ((grid.points %*% D) >= C) 227 wmask <- apply(wmask, 1, all) 228 wmask[wmask] <- 1 229 wmask[!wmask] <- NA 230 wmask <- matrix(wmask, ncol = ngrid) 231 } 232 233 w1 <- w1[rep(1:ngrid, ngrid), ] 234 w2 <- w2[rep(1:ngrid, each = ngrid), ] 235 ind.select <- function(i, w1, w2, selection) { 236 iset1 <- selection[i, ] 237 iset1 <- iset1[!is.na(iset1)] 238 iset2 <- !selection[i, ] 239 iset2 <- iset2[!is.na(iset2)] 240 (length(iset1) > 4) && (length(iset2) > 4) 241 sum((w1[i, iset1] > exp(-2)) & (w2[i, iset1] > exp(-2)), na.rm = TRUE) > 4 && 242 sum((w1[i, iset2] > exp(-2)) & (w2[i, iset2] > exp(-2)), na.rm = TRUE) > 4 243 } 244 245 beta1 <- beta1 * wmask 246 beta2 <- beta2 * wmask 247 var.angle <- atan2(beta2, beta1) + pi / 2 248 var.angle <- as.vector(var.angle) 249 250 sig2 <- sm.sigma(x, y) 251 A <- sig2$qmat 252 shat <- sig2$estimate 253 ts <- 0 254 B <- matrix(0, nrow = n, ncol = n) 255 256 for (ang in ((1:nangles) * pi / nangles)) { 257 258 angle <- var.angle 259 angle <- rep(ang, ngrid^2) * angle / angle 260 ev.points <- matrix(0, nrow=ngrid^2, ncol = 2) 261 ev.points[, 1] <- rep(eval.points[, 1], ngrid) 262 ev.points[, 2] <- rep(eval.points[, 2], rep(ngrid, ngrid)) 263 selection <- matrix(rep(cos(angle), n), ncol = n) * 264 (matrix(rep(x[, 2], ngrid^2), ncol = n, byrow = TRUE) 265 - matrix(rep(ev.points[, 2], n), ncol = n)) 266 selection <- selection - matrix(rep(sin(angle), n), ncol = n) * 267 (matrix(rep(x[, 1], ngrid^2), ncol = n, byrow = TRUE) 268 - matrix(rep(ev.points[, 1], n), ncol = n)) 269 270 selection <- (selection > 0) 271 ind <- sapply(1:nrow(selection), ind.select, 272 w1 = w1, w2 = w2, selection = selection) 273 ev.points <- ev.points[ind, ] 274 selection <- selection[ind, ] 275 276 selection[selection > 0] <- 1 277 selection[selection <= 0] <- -1 278 wl <- sm.discon.weight2(x, ev.points, h, (1 + selection) / 2) 279 wr <- sm.discon.weight2(x, ev.points, h, (1 - selection) / 2) 280 w <- wl - wr 281 if (doublesmooth) w <- sm.weight2(ev.points, ev.points, hd) %*% w 282 dhat <- as.vector(w %*% y) 283 s.e. <- as.vector(shat * sqrt((w^2) %*% rep(1,n))) 284 ts <- ts + sum((dhat / s.e.) ^2 ) 285 B <- B + t(w) %*% diag((shat/s.e.)^2) %*% w 286 } 287 288 C <- B - A * ts 289 p <- p.quad.moment(C, diag(n), 0, 0) 290 291 # Calculations for the reference band 292 293 angle <- var.angle 294 ev.points <- as.matrix(expand.grid(eval.points[, 1], eval.points[, 2])) 295 selection <- matrix(rep(cos(angle), n), ncol = n) * 296 (matrix(rep(x[, 2], ngrid^2), ncol = n, byrow = TRUE) 297 -matrix(rep(ev.points[, 2], n), ncol = n)) 298 selection <- selection - matrix(rep(sin(angle), n), ncol = n) * 299 (matrix(rep(x[, 1], ngrid^2), ncol = n, byrow = TRUE) 300 -matrix(rep(ev.points[, 1], n), ncol = n)) 301 selection <- (selection > 0) 302 ind <- sapply(1:nrow(selection), ind.select, 303 w1 = w1, w2 = w2, selection = selection) 304 ev.points <- ev.points[ind, ] 305 selection <- selection[ind, ] 306 307 selection[selection > 0] <- 1 308 selection[selection <= 0] <- -1 309 wl <- sm.discon.weight2(x, ev.points, h, (1 + selection) / 2) 310 wr <- sm.discon.weight2(x, ev.points, h, (1 - selection) / 2) 311 w <- wl - wr 312 if (doublesmooth) w <- sm.weight2(ev.points, ev.points, hd) %*% w 313 dhat <- as.vector(w %*% y) 314 s.e. <- as.vector(shat * sqrt((w^2) %*% rep(1, n))) 315 std <- rep(NA, ngrid * ngrid) 316 std[ind] <- dhat / s.e. 317 std <- matrix(abs(std), ncol = ngrid) 318 319 results <- list(p = p, sigma = shat, eval.points = eval.points, 320 st.diff = std, angle = matrix(angle, ncol = ngrid)) 321 322 if (opt$display != "none") { 323 if (!opt$add) 324 plot(x[, 1], x[, 2], xlab = opt$xlab, ylab = opt$ylab, 325 xlim = opt$xlim, ylim = opt$ylim, pch = opt$pch, col = opt$col.points) 326 mx <- max(std, na.rm = TRUE) 327 if (mx > 2.5) 328 contour(eval.points[, 1], eval.points[, 2], std, 329 levels = seq(2.5, mx, by = 0.5), 330 lty = opt$lty, col = opt$col, add = TRUE) 331 } 332 333 if (opt$verbose) 334 cat(paste("Test of continuity: significance = ", round(p, 3)), "\n") 335 336 invisible(results) 337 338 } 339 340#------------------------------------------------------------- 341 342sm.regression.invert <- function(a11,a12,a13,a22,a23,a33,c1,c2,c3) { 343 # Creates local linear intercept or slopes with two covariates 344 d <- a22 * a33 - a23^2 345 b1 <- 1 / (a11 - ((a12*a33 - a13*a23)*a12 + 346 (a13*a22 - a12*a23)*a13)/d) 347 b2 <- (a13*a23 - a12*a33) * b1 / d 348 b3 <- (a12*a23 - a13*a22) * b1 / d 349 est <- b1 * c1 + b2 * c2 + b3 * c3 350 invisible(est) 351 } 352 353#------------------------------------------------------- 354 355sm.discon.weight2 <- function(x, eval.points, h, selection, 356 weights = rep(1, nrow(x))) { 357 358# Amended version of sm.weight2 which uses different points 359# at each grid position 360 361 n <- nrow(x) 362 ne <- nrow(eval.points) 363 364 wd1 <- matrix(rep(eval.points[,1], rep(n, ne)), ncol = n, byrow = TRUE) 365 wd1 <- wd1 - matrix(rep(x[,1], ne), ncol = n, byrow = TRUE) 366 w <- exp(-.5 * (wd1 / h[1])^2) 367 wd2 <- matrix(rep(eval.points[,2], rep(n, ne)), ncol = n, byrow = TRUE) 368 wd2 <- wd2 - matrix(rep(x[,2], ne), ncol = n, byrow = TRUE) 369 w <- w * exp(-.5 * (wd2 / h[2])^2) 370 w <- w * matrix(rep(weights, ne), ncol = n, byrow = TRUE) 371 w <- w * selection 372 373 a11 <- w %*% rep(1, n) 374 a12 <- (w * wd1 ) %*% rep(1, n) 375 a13 <- (w * wd2 ) %*% rep(1, n) 376 a22 <- (w * wd1^2 ) %*% rep(1, n) 377 a23 <- (w * wd1 * wd2) %*% rep(1, n) 378 a33 <- (w * wd2^2 ) %*% rep(1, n) 379 380 d <- a22 * a33 - a23^2 381 382 b1 <- 1 / (a11 - ((a12*a33 - a13*a23)*a12 + (a13*a22 - a12*a23)*a13)/d) 383 b2 <- (a13*a23 - a12*a33) * b1 / d 384 b3 <- (a12*a23 - a13*a22) * b1 / d 385 386 wt <- matrix(rep(b1, n), ncol = n) 387 wt <- wt + matrix(rep(b2, n), ncol = n) * wd1 388 wt <- wt + matrix(rep(b3, n), ncol = n) * wd2 389 w <- wt * w 390 391 invisible(w) 392 393 } 394