# [R] S-code for piecewise regression

Date: Sat 05 Mar 2005 - 18:56:24 EST

"aicc" <- function(model)
{
# Calculates AICc for a nls model fit (e.g. pw.sharp)
# Reference: Anderson, D.R., K.P. Burnham and W.L. Thompson. 2002. Null
# hypothesis testing: problems, prevalence, and an alternative. Journal of Wildlife
# Management 64:912-923.
#
# model = nls model fit
#

```	sumfit <- summary(model)
N<-sum(sumfit\$df)
aicc <- N * log(sumfit\$df[2]*sumfit\$sigma^2/N) + 2 * sumfit\$df[1] + 2*
sumfit\$df[1]*( sumfit\$df[1]+1)/(N- sumfit\$df[1]-1)
return(aicc)
```

}

"aicc.fixedgamma" <- function(model)
{
# Calculates AICc for a nls model fit (e.g. pw.sharp) where gamma was fixed.
# Thus need to adjust for one extra estimted parameter (change the error and
# parameter df)
# Reference: Anderson, D.R., K.P. Burnham and W.L. Thompson. 2002. Null
# hypothesis testing: problems, prevalence, and an alternative. Journal of Wildlife
# Management 64:912-923.
#
# model = nls model fit
#

```	sumfit <- summary(model)
N <- sum(sumfit\$df)
K<-sumfit\$df[1]+1
aicc <- N * log((sumfit\$df[2] * sumfit\$sigma^2)/N) + 2 * K + (2 * K * (K +
1))/(N - K - 1)
return(aicc)
```

}
"basic.q" <- function(z, g)
{
# From Chiu 2002 Thesis
# one-parameter basic bent cable function

if(g > 0) ((z + g)^2/(4 * g)) * (z >= - g & z <= g) + z * (z > g) else z * (z > 0) }

"break.boot.ci" <- function(B, x, y, fnls, p, wts, ...) {
# Bootstrap the regression Data using the empirical distribution of the residuals.
# This function assumes that the variance of the errors is constant.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# B = number of "good" bootstrap samples required
# fnls = appropriate piecewise model function (e.g. pw.sharp)
# p = vector of starting value(s) required for the fnls function
# wts = vector of weights to use in fitting
#
# rs is the random seed. This can be used to replicate results.
# The original fit is in the first column of bootb
# Each row of failboot gives a bootstrap sample of y which did not
# converge to an answer when fitted with fnls. These models can
# be examined to diagnose problems in the bootstrapping.
#

```	rs <- .Random.seed
datadf <- data.frame(x = x, y = y, weights = wts)
e <- resid(data.nls)
bootb <- data.frame(coef(data.nls))
failboot <- NULL
yfit <- fitted(data.nls)
while(ncol(bootb) < (B + 1)) {
j <- sample(n, replace = T)	#resample residuals
d\$y <- (yfit + e[j])/sqrt(d\$weights)	#add to fitted y to give a bootstrapped y
boot.nls <- fnls(d\$x, d\$y, p, d\$weights, ...)
if(length(boot.nls) != 1)
bootb <- cbind(bootb, coef(boot.nls))
#if doesn't fail to converge, save coefficients
if(length(boot.nls) == 1)
failboot <- rbind(failboot, d\$y)	#if fails to converge, save bootstrapped y
}
return(rs, bootb, failboot)
```

}

"gen.lin.hyp" <- function(contrast, fit) {
# Runs a t-test for the hypothesis that a general linear combination
# of model parameters equals zero.
# contrast = vector representing the linear combination of parameters
# e.g. c(0,0,1,1) represents the combination 0*b0 + 0*b1 + 1*b2 + 1*b3 = 0
# fit = model fitted from a nls model (e.g. pw.sharp)
#

```	summ.fit <- summary(fit)
cov.mat <- summ.fit\$cov.unscaled * summ.fit\$sigma^2
se.comb <- sqrt(t(contrast) %*% cov.mat %*% contrast)
tc <- t(contrast) %*% fit\$param/se.comb
pvalue <- 2 * (1 - pt(abs(tc), summ.fit\$df[2]))
return(tc, pvalue)
```

}

"invert.ftest.bentcab" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a bent cable
# piecewise model.
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
# Gamma is unconstrained
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# p0 = vecot of initial estimate of breakpoint and gamma; e.g. p0=c(br0,g0)
# weights = vector of weights to use in fitting
# startstep = number of se(br) from brmax in which to search for the root above
# The nonlinear model is difficult to fit for fixed values of br far away from
# the max.
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	INF <- 1000000000
datadf\$dx - tau, gamm)), start = list(tau = p0[1], gamm = p0[2]), algorithm = "plinear")
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["tau"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["tau", "tau"] * ds2)
# start at brmax-startstep*se(br) and search towards brmax for root
n <- length(dy)
const <-  - error.df - qf(alpha, 1, error.df)
dxsort <- sort(unique(dx))
br <- brmax - startstep * sebr
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx, basic.q(datadf\$dx -
br, gamm)), start = list(gamm = p0[2]), algorithm = "plinear"))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brlower <- br
else {
br <- br + xeps
while((sign(obj1) == sign(obj2)) & (br < brmax)) {
obj1 <- obj2
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx,
basic.q(datadf\$dx - br, gamm)), start = list(gamm = p0[2]), algorithm =
"plinear"))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br + xeps
}
if(br >= brmax)
brlower <-  - INF
else brlower <- br - 2 * xeps
}
br <- brmax + startstep * sebr
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx, basic.q(datadf\$dx -
br, gamm)), start = list(gamm = p0[2]), algorithm = "plinear"))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brupper <- br
else {
br <- br - xeps
while((sign(obj1) == sign(obj2)) & (br > brmax)) {
obj1 <- obj2
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx,
basic.q(datadf\$dx - br, gamm)), start = list(gamm = p0[2]), algorithm =
"plinear"))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br - xeps
}
if(br <= brmax)
brupper <- INF
else brupper <- br + 2 * xeps
}
return(brmax, sebr, brlower, brupper)
```

}

"invert.ftest.bentcab.x" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95) {
#
# Special modification for the case where gamma is fixed.
#
# Inverted F-test confidence interval for the breakpoint in a bent cable
# piecewise model.
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# p0 = vecot of initial estimate of breakpoint and gamma; e.g. p0=c(br0,g0)
# weights = vector of weights to use in fitting
# startstep = number of se(br) from brmax in which to search for the root above
# The nonlinear model is difficult to fit for fixed values of br far away from
# the max.
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	INF <- 1000000000
datadf <- data.frame(dy, dx, weights, g = p0[2])
g <- p0[2]
nlsmax <- nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx, basic.q(dx - tau, g[1])), data =
datadf, start = list(tau = p[1]), algorithm = "plinear")
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["tau"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["tau", "tau"] * ds2)
# start at brmax-startstep*se(br) and search towards brmax for root
n <- length(dy)
const <-  - error.df - qf(alpha, 1, error.df)
dxsort <- sort(unique(dx))
br <- brmax - startstep * sebr
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brlower <- br
else {
br <- br + xeps
while((sign(obj1) == sign(obj2)) & (br < brmax)) {
obj1 <- obj2
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br + xeps
}
if(br >= brmax)
brlower <-  - INF
else brlower <- br - 2 * xeps
}
br <- brmax + startstep * sebr
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brupper <- br
else {
br <- br - xeps
while((sign(obj1) == sign(obj2)) & (br > brmax)) {
obj1 <- obj2
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br - xeps
}
if(br <= brmax)
brupper <- INF
else brupper <- br + 2 * xeps
}
return(brmax, sebr, brlower, brupper)
```

}

"invert.ftest.benthyp" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95) {
#
# find L.S. estimates of b0, b1, b2 and br in the bent hyperbola piecewise model
# dy = b0 + b1*(dx-br) + b2*sqrt((dx-br)^2+g^2/4)
# p0 = vector of initial estimates of br and g; e.g. p0=c(br0,g0)
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
# as alpha% confidence interval for br (i.e. alpha=.95 for 95% c.i.)
# Gamma is unconstrained
# S(beta(br), br) is the constrained SSerror when br has value br
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# startstep = number of se(br) from brmax in which to search for the root above
# The nonlinear model is difficult to fit for fixed values of br far away from
# the max.
# sample call:
#

```	INF <- 1000000000
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
# start at brmax-startstep*se(br) and search towards brmax for root
n <- length(dy)
const <-  - error.df - qf(alpha, 1, error.df)
dxsort <- sort(unique(dx))
br <- brmax - startstep * sebr
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br, sqrt((distance -
br)^2 + g^2/4)), start = list(g = p0[2]), algorithm = "plinear"))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brlower <- br
else {
br <- br + xeps
while((sign(obj1) == sign(obj2)) & (br < brmax)) {
obj1 <- obj2
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br,
sqrt((distance - br)^2 + g^2/4)), start = list(g = p0[2]), algorithm =
"plinear"))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br + xeps
}
if(br >= brmax)
brlower <-  - INF
else brlower <- br - 2 * xeps
}
br <- brmax + startstep * sebr
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br, sqrt((distance -
br)^2 + g^2/4)), start = list(g = p0[2]), algorithm = "plinear"))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brupper <- br
else {
br <- br - xeps
while((sign(obj1) == sign(obj2)) & (br > brmax)) {
obj1 <- obj2
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br,
sqrt((distance - br)^2 + g^2/4)), start = list(g = p0[2]), algorithm =
"plinear"))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br - xeps
}
if(br <= brmax)
brupper <- INF
else brupper <- br + 2 * xeps
}
return(brmax, sebr, brlower, brupper)
```

}

"invert.ftest.sharp" <- function(dx, dy, br0, weights, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a sharp piecewise model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when br has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	INF <- 1000000000
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
# need to start at the 2nd smallest observed dx and search for root
n <- length(dy)
const <-  - error.df - qf(alpha, 1, error.df)
dxsort <- sort(unique(dx))
br <- dxsort[2]
sum.fit <- summary(lm(datadf\$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights = weights))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brlower <- br
else {
br <- br + xeps
while((sign(obj1) == sign(obj2)) & (br < brmax)) {
obj1 <- obj2
sum.fit <- summary(lm(datadf\$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights
= weights))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br + xeps
}
if(br >= brmax)
brlower <-  - INF
else brlower <- br - 2 * xeps
}
br <- rev(dxsort)[2]
sum.fit <- summary(lm(datadf\$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights = weights))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brupper <- br
else {
br <- br - xeps
while((sign(obj1) == sign(obj2)) & (br > brmax)) {
obj1 <- obj2
sum.fit <- summary(lm(datadf\$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights
= weights))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br - xeps
}
if(br <= brmax)
brupper <- INF
else brupper <- br + 2 * xeps
}
return(brmax, sebr, brlower, brupper)
```

}

"invert.ftest.tanh" <- function(dx, dy, p0, weights, startstep = 4, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a hyperbolic tangent
# piecewise model.
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <= F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
# Gamma is unconstrained
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# p0 = vecot of initial estimate of breakpoint and gamma; e.g. p0=c(br0,g0)
# weights = vector of weights to use in fitting
# startstep = number of se(br) from brmax in which to search for the root above
# The nonlinear model is difficult to fit for fixed values of br far away from
# the max.
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	INF <- 1000000000
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
# start at brmax-startstep*se(br) and search towards brmax for root
n <- length(dy)
const <-  - error.df - qf(alpha, 1, error.df)
dxsort <- sort(unique(dx))
br <- brmax - startstep * sebr
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br, (dx - br) *
tanh((dx - br)/g)), start = list(g = p0[2]), algorithm = "plinear"))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brlower <- br
else {
br <- br + xeps
while((sign(obj1) == sign(obj2)) & (br < brmax)) {
obj1 <- obj2
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br,
(dx - br) * tanh((dx - br)/g)), start = list(g = p0[2]), algorithm =
"plinear"))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br + xeps
}
if(br >= brmax)
brlower <-  - INF
else brlower <- br - 2 * xeps
}
br <- brmax + startstep * sebr
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br, (dx - br) *
tanh((dx - br)/g)), start = list(g = p0[2]), algorithm = "plinear"))
obj1 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
obj2 <- obj1
if(obj1 == 0)
brupper <- br
else {
br <- br - xeps
while((sign(obj1) == sign(obj2)) & (br > brmax)) {
obj1 <- obj2
sum.fit <- summary(nls(sqrt(weights) * datadf\$dy ~ sqrt(weights) * cbind(1, dx - br,
(dx - br) * tanh((dx - br)/g)), start = list(g = p0[2]), algorithm =
"plinear"))
obj2 <- (sum.fit\$sigma^2 * sum.fit\$df[2])/ds2 + const
br <- br - xeps
}
if(br <= brmax)
brupper <- INF
else brupper <- br + 2 * xeps
}
return(brmax, sebr, brlower, brupper)
```

}

"lack.of.fit" <- function(x, y, fit, weights) {
# tests of lack-of-fit for nonlinear models versus saturated model
#
# x = vector of x data to fit
# y = vector of y data to fit
# fit = model fit to be tested fot lack of fit (e.g. resulting from pw.sharp)
# weights = vector of weights to use in fitting
#

```	factor.fit <- lm(y ~ factor(x), weights = weights)
dfpe <- factor.fit\$df.residual
sspe <- (summary(factor.fit))\$sigma^2 * dfpe
sse <- (summary(fit))\$sigma^2 * (summary(fit))\$df[2]
dflof <- (summary(fit))\$df[2] - dfpe
Fc <- ((sse - sspe) * dfpe)/dflof/sspe
pvalue <- 1 - pf(Fc, dflof, dfpe)
return(Fc, pvalue)
```

}

"llsurface.bentcab" <- function(dx, dy, br, g, startvals, weights) {
# Calculates a linear fit for each grid point (grid of possible br values).
# Calculates a log-likelihood surface across the grid.
# For the bent cable piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
# persp(out\$br,out\$gamma,out\$loglik,xlab="Value of the breakpoint, br",
# ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting values for breakpoint and gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting
#

```	y <- dy
x <- dx
datadf <- data.frame(x = dx, y = dy, weights = weights)
p0 <- startvals
nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, distance, basic.q(distance - br, g)),
start = list(br = p0[1], g = p0[2]), data = datadf, algorithm = "plinear")
nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]	#sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)

#n <- length(y)	#const <- 4. - n - qf(alpha, 1., n - 4.)

out <- NA
out\$loglik <- matrix(nrow = length(br), ncol = length(g))
# walk through all values of br and gamma, calculating the deviance
for(i in 1:length(br)) {
for(j in 1:length(g)) {
out\$loglik[i, j] <- summary(lm(datadf\$y ~ cbind(distance, basic.q(distance - br[i],
g[j])), weights = weights))\$sigma^2
}
}
out\$loglik <- (out\$loglik * (error.df + 2) - ds2 * error.df)/ds2/2
out\$br <- br
out\$gamma <- g
return(out)
```

}

"llsurface.bentcab.x" <- function(dx, dy, br, g, startvals, weights) {
#
# Special modification for the case where gamma is fixed.
#
# Calculates a linear fit for each grid point (grid of possible br values).
# Calculates a log-likelihood surface across the grid.
# For the bent cable piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
# persp(out\$br,out\$gamma,out\$loglik,xlab="Value of the breakpoint, br",
# ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting value for breakpoint, fixed value for gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting
#

```	y <- dy
x <- dx
p0 <- startvals
g0 <- p0[2]	#datadf <- data.frame(x = dx, y = dy, weights = weights, p0 = startvals)
datadf <- data.frame(x = dx, y = dy, weights = weights, g0)
#nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1., x, basic.q(x - br, p0[2.])), start #= list(br = p0[1.]), data = datadf, algorithm = "plinear")
nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, basic.q(x - br, g0[1])), start = list(
br = p0[1]), data = datadf, algorithm = "plinear")
nlsum <- summary(nlsmax)
ds2 <- nlsum\$sigma^2	#sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)

#n <- length(y) #const <- 4. - n - qf(alpha, 1., n - 4.)

error.df <- nlsum\$df[2]
out <- NA	# walk through all values of br and gamma, calculating the deviance
out\$loglik <- matrix(nrow = length(br), ncol = length(g))
for(i in 1:length(br)) {
for(j in 1:length(g)) {
g[j])), weights = weights))\$sigma^2
}
}
out\$loglik <- (out\$loglik * (error.df + 1) - ds2 * error.df)/ds2/2
out\$br <- br
out\$gamma <- g
return(out)
```

}

"llsurface.benthyp" <- function(dx, dy, br, g, startvals, weights) {
# Calculates a linear fit for each grid point (grid of possible br values).
# Plots a log-likelihood surface across the grid.
# For the bent hyperbola piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
# persp(out\$br,out\$gamma,out\$loglik,xlab="Value of the breakpoint, br",
# ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting values for breakpoint and gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting
#

```	y <- dy
x <- dx
datadf <- data.frame(x = dx, y = dy, weights = weights)
p0 <- startvals
nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), sqrt((distance - br)^2 +
g^2/4)), start = list(br = p0[1], g = p0[2]), data = datadf, algorithm = "plinear")
nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]	#sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)

#n <- length(y)	#const <- 4. - n - qf(alpha, 1., n - 4.)

out <- NA
out\$loglik <- matrix(nrow = length(br), ncol = length(g))
# walk through all values of br and gamma, calculating the deviance
for(i in 1:length(br)) {
for(j in 1:length(g)) {
out\$loglik[i, j] <- summary(lm(datadf\$y ~ cbind((distance - br[i]), sqrt((distance -
br[i])^2 + g^2/4)), weights = weights))\$sigma^2
}
}
out\$loglik <- (out\$loglik * (error.df + 2) - ds2 * error.df)/ds2/2
out\$br <- br
out\$gamma <- g
return(out)
```

}

"llsurface.sharp" <- function(dx, dy, br, br.start, weights) {
# Calculates a linear fit for each grid point (grid of possible br values).
# Calculates a log-likelihood surface across the grid.
# For the sharp piecewise regression model.
#
# Can use to create a perspective log-likelihood plot:
# persp(out\$br,out\$gamma,out\$loglik,xlab="Value of the breakpoint, br",
# ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# br.start = initial starting values for breakpoint; e.g. br.start=br0
# weights = vector of weights to use in fitting
#

```	y <- dy
x <- dx
datadf <- data.frame(x = dx, y = dy, weights = weights)
br0 <- br.start
nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, ifelse(x > br, x - br, 0)), start =
list(br = br0), control = list(maxiter = 15, tolerance = 1e-005, minscale = 1e-005), data
nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]	#sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)

#n <- length(y) #const <- 4. - n - qf(alpha, 1., n - 4.)

out <- data.frame(br = br, loglik = rep(NA, length(br)))
# walk through all values of br, calculating the deviance
out\$loglik <- sapply(br, loopfn.sharp, datadf, weights)
out\$loglik <- (out\$loglik * (error.df + 1) - ds2 * error.df)/ds2/2
print(cbind(br.fitted = brmax, sebr.fitted = sebr, sigmasq.fitted = ds2))
return(out)
```

}

"llsurface.tanh" <- function(dx, dy, br, g, startvals, weights) {
# Calculates a linear fit for each grid point (grid of possible br, gamma)
# Plots a scaled deviance surface across the grid,
# for the hyperbolic tangent piecewise regression model.
# This surface can be used to compute approx. conf regions by comparing with
# quantiles of the F(2, nlsum\$df[2]), since the scaling incorporates the 2 df
# for the 2 parameters br and gamma.
#
# Can use to create a perspective log-likelihood plot:
# persp(out\$br,out\$gamma,out\$loglik,xlab="Value of the breakpoint, br",
# ylab="Value of gamma, g", zlab="Scaled log-likelihood")
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br = vector of grid points of the breakpoint parameter to be used
# g = vector of grid points of gamma to be used
# startvals = initial starting values for breakpoint and gamma; e.g. startvals=c(br0,g0)
# weights = vector of weights to use in fitting; proportional to 1/Var(dy)
#

```	y <- dy
x <- dx
datadf <- data.frame(x = dx, y = dy, weights = weights)
p0 <- startvals
nlsmax <- nls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), (distance - br) * tanh((
distance - br)/g)), start = list(br = p0[1], g = p0[2]), data = datadf, algorithm =
"plinear")
nlsum <- summary(nlsmax)	#brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]	#sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)

#n <- length(y)	#const <- 4. - n - qf(alpha, 1., n - 4.)

out <- NA
out\$loglik <- matrix(nrow = length(br), ncol = length(g))
# walk through all values of br and gamma, calculating the deviance
for(i in 1:length(br)) {
for(j in 1:length(g)) {
out\$loglik[i, j] <- summary(lm(datadf\$y ~ cbind((distance - br[i]), (distance - br[
i]) * tanh((distance - br[i])/g[j])), weights = weights))\$sigma^2
}
}
out\$loglik <- (out\$loglik * (error.df + 2) - ds2 * error.df)/ds2/2
out\$br <- br
out\$gamma <- g
return(out)
```

}

"loopfn.sharp" <- function(br, data, weights) {
# internal function for llsurface.sharp
#

```	datadf <- data
br0 <- br
loglik <- summary(lm(datadf\$y ~ cbind(x, ifelse(x > br0, x - br0, 0)), weights = weights))\$sigma^2
return(loglik)
```

}

"mynls" <- function(formula, data = sys.parent(), start = if(!is.null(attr(data, "parameters"))) attr(data,

"parameters") else nls.initial(), control, algorithm = "default", trace = F) {

```	convert.twiddle <- function(formula)
{
if(length(formula) < 3)
return(formula)
form <- call("~", call("-", formula[[2]], formula[[3]]))
attr(form, "class") <- "formula"
form
}
if(is.numeric(data))
data <- sys.frame(data)
cl <- class(data)
if(inherits(data, "pframe")) {
class(data) <- NULL
pp <- parameters(data)
if(length(pp)) {
np <- names(pp)
if(any(match(np, names(data), 0)))
stop("can't have variables, parameters with\nsame na\nme")
data[np] <- pp
}
}
else if(inherits(data, "data.frame"))
cl <- c("pframe", cl)
class(data) <- NULL	## First, figure out the par. names, make start a list
switch(mode(start),
numeric = {
.parameters <- names(start)
start <- as.list(start)
names(start) <- .parameters
}
,
list = {
.parameters <- names(start)
}
,
NULL = .parameters <- parameter.names(formula, data),
stop("\"start\" should be numeric or list"))
if(!length(.parameters))
stop("names for parameters needed, from formula or from start")
pn <- .parameters
.expr <- formula	## select the algorithm and possibly massage the formula
if(length(start)) data[.parameters] <- start	#used by setup_nonlin
data\$.parameters <- .parameters
nl.frame <- new.frame(data, F)
frame.attr("class", nl.frame) <- cl	# in case data is returned
formula <- switch(algorithm,
plinear = {
if(length(formula) < 3)
stop("formula for plinear algorithm be of the form\nresp ~ array")
response <- eval(formula[[2]], nl.frame)
design <- eval(formula[[3]], nl.frame)
nnobs <- length(response)
nlinear <- if(is.matrix(design)) dim(design)[2] else length(design)/nnobs
form <- call("~", formula[[3]])
attr(form, "class") <- "formula"
form
}
,
default = convert.twiddle(formula))
dims <- .C("setup_nonlin",
n = integer(3),
list(formula),
as.integer(nl.frame))\$n
npar <- dims[1]
nderiv <- dims[2]
nobs <- dims[3]
resp <- eval(.expr[[2]], nl.frame)
if(length(.expr) == 2) resp[] <- 0	## setup_nonlin will have set .parameters if missing
if(is.null(start)) {
start <- list()
if(is.null(names(start)) && length(start) == length(pn))
names(start) <- pn
for(i in .parameters)
start[[i]] <- get(i, frame = nl.frame)
}
asgn <- start
si <- 1
for(i in 1:length(asgn)) {
ni <- length(asgn[[i]])
asgn[[i]] <- seq(from = si, length = ni)
si <- si + ni
}
start <- unlist(start)
controlvals <- nls.control()
if(!missing(control))
controlvals[names(control)] <- control
ret.data <- is.logical(ret.data <- controlvals\$data) && ret.data
max.iterations <- if(is.null(controlvals\$maxiter)) 50 * npar else controlvals\$maxiter
settings <- c(0, max.iterations, controlvals\$minscale, controlvals\$tolerance, 0)
if(algorithm == "plinear") {
settings[1] <- 1
dims <- c(dims, nlinear)
start <- c(start, numeric(nlinear))
pn <- c(pn, paste(".lin", 1:nlinear, sep = ""))
dims[3] <- nnobs
outmat <- array(c(response, numeric(nnobs * (npar + nlinear))), c(nnobs, npar + nlinear + 1
))
}
else outmat <- array(0, c(nobs, npar + 1))
storage.mode(outmat) <- "double"
storage.mode(start) <- "double"
nls.trace <- if(missing(trace)) controlvals\$trace else trace
std.trace <- FALSE
if(is.logical(nls.trace)) {
if(std.trace <- nls.trace)
nls.trace <- "trace.nls"
else nls.trace <- NULL
}
else std.trace <- is.character(nls.trace) && nls.trace == "trace.nls"
if(std.trace) {
assign("trace.mat", array(0, c(max.iterations, npar + 2), list(NULL, c("obj.", "conv.",
paste("par", 1:npar)))), frame = 1)
assign("trace.expr", expression(trace.mat[last.iteration,  ] <- it.row), frame = 1)
}
z <- .C("do_nls",
parameters = start,
dims = as.integer(dims),
control = as.double(settings),
outmat = outmat,
trace = list(nls.trace))
if(z\$control[5]) nls.out <- c("skip") else {
nls.out <- list(parameters = z\$parameters, formula = .expr, call = match.call(), residuals
= z\$outmat[, 1])
if(algorithm == "plinear") {
class(nls.out) <- c("nls.pl", "nls")
R <- qr(z\$outmat[, -1])\$qr[1:(npar + nlinear),  , drop = F]
}
else {
class(nls.out) <- "nls"
R <- qr(z\$outmat[, -1])\$qr[1:npar,  , drop = F]
}
R[lower.tri(R)] <- 0
nls.out\$R <- R
nls.out\$fitted.values <- resp - nls.out\$residuals
nls.out\$assign <- asgn
if(ret.data) {
data <- sys.frame(nl.frame)
data\$.parameters <- NULL
#dbdetach does the right thing--only matters that nl.frame>1
if(inherits(data, "pframe"))
data <- dbdetach(data, nl.frame)
nls.out\$data <- data
}
if(std.trace && exists("last.iteration", frame = 1))
nls.out\$trace <- get("trace.mat", frame = 1)[1:get("last.iteration", frame = 1),  ]
}

#       stop(switch(as.integer(z\$control[5]),
#               "step factor reduced below minimum",
#               "maximum number of iterations exceeded",
#               "singular design matrix",

nls.out
```

}

"plot.invert.ftest.bentcab" <- function(dx, dy, br0, g0, weights, limits, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a bent cable piecewise
# model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# limits = vector length 2 of start and end limits
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	datadf <- data.frame(dy, dx, weights, br = br0)
p0 <- c(br0, g0)
print(nlsmax)
nlsum <- summary(nlsmax)
print(nlsum)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
print(brmax)
Fcrit <- qf(alpha, 1, error.df)	# undx <- sort(unique(dx))

# brgrid <- seq(undx[2], rev(undx)[2], xeps)

brgrid <- seq(limits[1], limits[2], xeps)
obj1 <- NA
for(br in brgrid) {
print(br)
sum.fit <- summary(nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx, basic.q(dx - br,
g0)), start = list(g0 = p0[2]), data = datadf, algorithm = "plinear"))
obj1 <- cbind(obj1, sum.fit\$sigma)
}
Fstats <- (obj1[-1]^2 * sum.fit\$df[2])/ds2 - error.df
Fstats <- cbind(brgrid, Fstats)
plot(Fstats[, 1], Fstats[, 2], type = "l")
abline(h = Fcrit)
return(brmax, sebr, Fcrit, Fstats)
```

}

"plot.invert.ftest.bentcab.x" <- function(dx, dy, br0, g0, weights, xeps = 2, alpha = 0.95) {
# Special modification for the case where gamma is fixed.
#
# Inverted F-test confidence interval for the breakpoint in a bent cable piecewise
# model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	p0 <- c(br0, g0)
datadf <- data.frame(dy, dx, weights, p0)
nlsmax <- nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx, basic.q(dx - br, p0[2])), data =
datadf, start = list(br = p0[1]), algorithm = "plinear")
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
Fcrit <- qf(alpha, 1, error.df)
undx <- sort(unique(dx))
brgrid <- seq(undx[2], rev(undx)[2], xeps)
obj1 <- NA
for(br in brgrid) {
sum.fit <- summary(lm(datadf\$dy ~ cbind(distance, basic.q(distance - br, p0[2])), weights
= weights))
obj1 <- cbind(obj1, sum.fit\$sigma)
}
Fstats <- (obj1[-1]^2 * sum.fit\$df[2])/ds2 - error.df
Fstats <- cbind(brgrid, Fstats)
plot(Fstats[, 1], Fstats[, 2], type = "l")
abline(h = Fcrit)
return(brmax, sebr, Fcrit, Fstats)
```

}

"plot.invert.ftest.benthyp" <- function(dx, dy, br0, g0, weights, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a bent hyperbola piecewise
# model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	datadf <- data.frame(dy, dx, weights)
p0 <- c(br0, g0)
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
Fcrit <- qf(alpha, 1, error.df)
undx <- sort(unique(dx))
brgrid <- seq(undx[2], rev(undx)[2], xeps)
obj1 <- NA
for(br in brgrid) {
sum.fit <- summary(lm(datadf\$dy ~ cbind(distance - br, sqrt((distance - br)^2 + (g0)^2/4)),
weights = weights))
obj1 <- cbind(obj1, sum.fit\$sigma)
}
Fstats <- (obj1[-1]^2 * sum.fit\$df[2])/ds2 - error.df
Fstats <- cbind(brgrid, Fstats)
plot(Fstats[, 1], Fstats[, 2], type = "l")
abline(h = Fcrit)
return(brmax, sebr, Fcrit, Fstats)
```

}

"plot.invert.ftest.sharp" <- function(dx, dy, br0, weights, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a sharp piecewise
# model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2 <=F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# weights = vector of weights to use in fitting
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	datadf <- data.frame(dy, dx, weights)
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
Fcrit <- qf(alpha, 1, error.df)
undx <- sort(unique(dx))
brgrid <- seq(undx[2], rev(undx)[2], xeps)
obj1 <- NA
for(br in brgrid) {
sum.fit <- summary(lm(datadf\$dy ~ cbind(dx, ifelse(dx > br, dx - br, 0)), weights = weights
))
obj1 <- cbind(obj1, sum.fit\$sigma)
}
Fstats <- (obj1[-1]^2 * sum.fit\$df[2])/ds2 - error.df
Fstats <- cbind(brgrid, Fstats)
plot(Fstats[, 1], Fstats[, 2], type = "l")
abline(h = Fcrit)
return(brmax, sebr, Fcrit, Fstats)
```

}

"plot.invert.ftest.tanh" <- function(dx, dy, br0, g0, weights, limits, xeps = 2, alpha = 0.95) {
# Inverted F-test confidence interval for the breakpoint in a hyperbolic
# tangent piecewise
# model
# Find roots of [S(beta(br),br) - S(beta_hat, br_hat)]/ds2
# <=F(alpha,1,error.df)
# as alpha% confidence interval for br
# S(beta(br), br) is the constrained SSerror when breakpoint has value br
#
# dx = vector of the x data to fit
# dy = vector of the y data to fit
# br0 = initial estimate of breakpoint
# g0 = initial estimate of gamma
# weights = vector of weights to use in fitting
# limits = vector length 2 of start and end limits
# xeps = step size in grid search, smaller values are more precise but
# computationally intensive
# alpha=.95 for 95% confidence interval
#
# brmax = maximum likelihood (least squares) estimate of the breakpoint
# sebr = standard error of brmax
# brlower = lower estimate of confidence interval
# brupper = upper estimate of confidence interval
#

```	datadf <- data.frame(dy, dx, weights, br = br0)
p0 <- c(br0, g0)
nlsum <- summary(nlsmax)
brmax <- (coef(nlsmax))["br"]
ds2 <- nlsum\$sigma^2
error.df <- nlsum\$df[2]
sebr <- sqrt(nlsum\$cov.unscaled["br", "br"] * ds2)
print(brmax)
Fcrit <- qf(alpha, 1, error.df)	#       undx <- sort(unique(dx))

#       brgrid <- seq(undx[2], rev(undx)[2], xeps)

brgrid <- seq(limits[1], limits[2], xeps)
obj1 <- NA
for(br in brgrid) {
print(br)
sum.fit <- summary(nls(sqrt(weights) * dy ~ sqrt(weights) * cbind(1, dx - br, (dx - br) *
tanh((dx - br)/g0)), start = list(g0 = p0[2]), data = datadf, algorithm = "plinear"
))
obj1 <- cbind(obj1, sum.fit\$sigma)
}
Fstats <- (obj1[-1]^2 * sum.fit\$df[2])/ds2 - error.df
Fstats <- cbind(brgrid, Fstats)
plot(Fstats[, 1], Fstats[, 2], type = "l")
abline(h = Fcrit)
return(brmax, sebr, Fcrit, Fstats)
```

}

"pw.bentcab" <- function(x, y, p, weights, maxit = 50) {
# Fits bent cable piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting values (breakpoint, gamma)
# weights = vector of weights to use in fitting
#

```	datadf <- data.frame(x = x, y = y, weights = weights)
data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, basic.q(x - br, g)), data =
datadf, start = list(br = p[1], g = p[2]), algorithm = "plinear", control = list(maxiter =
maxit))
return(data.nls)
```

}

"pw.bentcab.x" <- function(x, y, p, weights, maxit = 50) {
# Fits bent cable piecewise regression with one breakpoint.
# Modification for the case when gamma needs to be fixed.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting value of breakpoint, fixed value for gamma
# weights = vector of weights to use in fitting
#

```	datadf <- data.frame(x = x, y = y, weights = weights, g = p[2])
data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, x, basic.q(x - br, g[1])), data =
datadf, start = list(br = p[1]), algorithm = "plinear", control = list(maxiter = maxit))
return(data.nls)
```

}

"pw.benthyp" <- function(x, y, p, weights, maxit = 50) {
# Fits bent hyperbola piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting values (breakpoint, gamma)
# weights = vector of weights to use in fitting
#

```	datadf <- data.frame(x = x, y = y, weights = weights)
data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), sqrt((distance - br
)^2 + g^2/4)), data = datadf, start = list(br = p[1], g = p[2]), algorithm = "plinear",
control = list(maxiter = maxit))
return(data.nls)
```

}

"pw.sharp" <- function(x, y, p, weights, minscale = 0.001, toler = 0.001, trace = F) {
# Fits sharp piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = starting value of breakpoint
# weights = vector of weights to use in fitting
#

```	datadf <- data.frame(x = x, y = y, weights = weights)
data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, distance, ifelse(distance > br,
distance - br, 0)), data = datadf, start = list(br = p), algorithm = "plinear", control = c(
minscale = minscale, tolerance = toler), trace = trace)
return(data.nls)
```

}

"pw.tanh" <- function(x, y, p, weights, maxit = 50) {
# Fits hyperbolic tangent piecewise regression with one breakpoint.
#
# x = vector of the x data to fit
# y = vector of the y data to fit
# p = vector of starting values (breakpoint, gamma)
# weights = vector of weights to use in fitting
#

```	datadf <- data.frame(x = x, y = y, weights = weights)
data.nls <- mynls(sqrt(weights) * y ~ sqrt(weights) * cbind(1, (distance - br), (distance - br) *
tanh((distance - br)/g)), data = datadf, start = list(br = p[1], g = p[2]), algorithm =
"plinear", control = list(maxiter = maxit))
return(data.nls)
```

}

R-help@stat.math.ethz.ch mailing list