I am interested in estimating a break point in response to one explanatory variable, while also including other variables as terms in the linear model. One reason to include additional terms is to account for spatial dependency/autocorrelation. The following code was adapted from code provided by @jbowman here, and also referred to here. It finds the break point for the model lm(Rare25~epi.i). What I would like to do is find a break point for Rare25 at a value of epi.i, while including the variable l.dist.mrc in a model such as lm(Rare25 ~ epi.i + l.dist.mrc).
Data:
epi.i<-c(7, 8, 9, 10, 8.01, 6, 7.01, 5.999, 4.003, 12, 2.01, 7.02, 10.01, 8.2, 5.9, 3.9, 6.999, 4.0001, 3.99, 6.001, 8.001, 5.99, 7.9, 6.99, 3.98)
Rare25<-c(4.471429, 4.551474, 4.204894, 4.456710, 3.504348, 4.175824, 4.298193, 4.406838, 3.058707, 4.451128, 1.000000, 4.327893, 3.580541, 4.082432, 3.630869, 4.352130, 3.919075, 3.795205, 3.380952, 2.993347, 3.775886, 3.766723, 3.852396, 3.923977, 4.308840)
l.dist.mrc<-c(2.7839, 2.7839, 2.7839, 3.0232, 3.0232, 2.5051, 2.5051, 2.7041, 2.9965, 3.0773, 3.4457, 3.1202, 3.1936, 3.4362, 3.3718, 2.7641, 2.7641, 3.7904, 3.4457, 3.7904, 3.1202, 2.5302, 2.5302 ,3.1242, 3.1241)
Code:
#Basis functions
bp = 4
b1 <- function(x, bp) ifelse(x < bp, bp - x, 0)
b2 <- function(x, bp) ifelse(x < bp, 0, x - bp)
#Wrapper
foo <- function(bp)
{
mod <- lm(Rare25 ~ b1(epi.i, bp) + b2(epi.i, bp))
deviance(mod)
}
search.range <- c(min(epi.i),max(epi.i))
foo.opt <- optimize(foo, interval = search.range)
bp <- foo.opt$minimum
bp
# confidence interval
foo.root <- function(bp, tgt)
{
foo(bp) - tgt
}
tgt <- foo.opt$objective + qchisq(0.95,1)
lb95 <- uniroot(foo.root, lower=(search.range[1]), upper=bp, tgt=tgt)
ub95 <- uniroot(foo.root, lower=bp, upper=search.range[2], tgt=tgt)
lb95$root
ub95$root
