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