[R] best way to fit a model

[Achim Zeileis]

On Fri, 9 Sep 2005 14:55:04 -0300 Ronaldo Reis-Jr. wrote:

> Hi,
> 
> I have some data that have this behaviour:
> 
> |
> |*******
> |        *
> |          * 
> |            *
> |              *
> |----------------
> 
> What is the best and simpler way to fit this in R?

If the changepoint is known, then this is straightforward using lm:

## generate example data
set.seed(20050909)
x <- seq(0, 10, by = 0.25)
y.mean <- ifelse(x >= 5, x, 5)
y <- y.mean + rnorm(41)
plot(y ~ x)
lines(y.mean ~ x)

## fit linear model with break
fm <- lm(y ~ I((x-5) * (x >= 5)))
fm
y.fit1 <- fitted(fm)
lines(y.fit1 ~ x, col = 2)

If it is unknown, it can be estimated using Vito Muggeo's segmented
package:

## estimate change point in x
library("segmented")
sfm <- segmented(lm(y ~ x), x, 6)
sfm
y.fit2 <- fitted(sfm)
lines(y.fit2 ~ x, col = 3)

This fits a continuous mean function. Alternatively, breakpoints() in
strucchange can be used to estimate a break point:

## estimate break point in x
library("strucchange")
bp <- breakpoints(y ~ x)
summary(bp)
y.fit3 <- fitted(bp)
lines(y.fit3 ~ x, col = 4)

This does not enforce that the line is continuous, hence the jump in the
fitted mean. Of course, the estimated breakpoint could be used to fit a
continuous line model, but this is not what is optimized in
breakpoints().
Z

> Thanks
> Ronaldo
> -- 
> Ela pilotava um Continenal 2001 com ignição automática Magiclic...
> --
> |>   // | \\   [***********************************]
> |   ( õ   õ )  [Ronaldo Reis Júnior                ]
> |>      V      [UFV/DBA-Entomologia                ]
> |    /     \   [36570-000 Viçosa - MG              ]
> |>  /(.''`.)\  [Fone: 31-3899-4007                 ]
> |  /(: :'  :)\ [chrysopa at insecta.ufv.br            ]
> |>/ (`. `'` ) \[ICQ#: 5692561 | LinuxUser#: 205366 ]
> |    ( `-  )   [***********************************]
> |>>  _/   \_Powered by GNU/Debian Woody/Sarge
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>