Dear all
I am trying to find one breakpoint in my data. I use segmented package.
Here is my data> dput(temp)
temp <- structure(list(spotreba = 0:40, sqsp = c(0, 1, 1.4142135623731,
1.73205080756888, 2, 2.23606797749979, 2.44948974278318, 2.64575131106459,
2.82842712474619, 3, 3.16227766016838, 3.3166247903554, 3.46410161513775,
3.60555127546399, 3.74165738677394, 3.87298334620742, 4, 4.12310562561766,
4.24264068711928, 4.35889894354067, 4.47213595499958, 4.58257569495584,
4.69041575982343, 4.79583152331272, 4.89897948556636, 5, 5.09901951359278,
5.19615242270663, 5.29150262212918, 5.3851648071345, 5.47722557505166,
5.56776436283002, 5.65685424949238, 5.74456264653803, 5.8309518948453,
5.91607978309962, 6, 6.08276253029822, 6.16441400296898, 6.2449979983984,
6.32455532033676), variable = structure(c(5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L,
5L, 5L, 5L, 5L), .Label = c("vod1", "vod2",
"vod3", "vod4", "vod5"
), class = "factor"), value = c(420, 438, 449, 456, 466, 472,
476, 481, 487, 496, 507, 520, 535, 554, 577, 606, 639, 677, 722,
776, 836, 897, 959, 1030, 1105, 1182, 1256, 1346, 1432, 1519,
1613, 1704, 1802, 1894, 1998, 2097, 2189, 2290, 2393, 2502, 2608
), fit = c(334.519176009225, 405.569071381122, 434.998901649351,
457.581204665804, 476.618966753019, 493.391571855035, 508.555165948505,
522.499529840424, 535.478627289478, 547.668862124917, 559.198672901077,
570.165020351817, 580.643233322383, 590.693216888952, 600.363541857,
609.694237534363, 618.718757496814, 627.465399316641, 635.958352929605,
644.218489884464, 704.452390602657, 811.399180061218, 915.82851606979,
1017.91022567874, 1117.79585929747, 1215.6213312916, 1311.50908941826,
1405.56991195017, 1497.904407726, 1588.60427702993, 1677.75337832018,
1765.42863614133, 1851.70081819772, 1936.63520392068, 2020.29216249248,
2102.72765487767, 2183.99367172737, 2264.13861689264, 2343.20764458285,
2421.2429568385, 2498.28406688275)), row.names = 165:205, class
"data.frame")
fit <- lm(value~sqsp, temp)
fit.s <- segmented(fit, seg.Z= ~sqsp, npsi=1, psi = 2, nboot=50, tol=1e-10,
quant=F)
plot(temp$sqsp, temp$value)
plot(fit.s, add=T)
You can see that first part is rather off measured points. I tried to fiddle
with seg.control but was not able to achieve better result. The only way how
to make closer fit to first part of data which I was able to find was to get
rid of points near estimated psi and refit model.
temp.v <- temp[-c(15:25),]
fit <- lm(value~sqsp, temp.v)
fit.s <- segmented(fit, seg.Z= ~sqsp, npsi=1, psi = 2, nboot=50, tol=1e-10,
quant=F)
plot(fit.s, add=T, col=2)
Is there any more elegant way how to achieve similar result?
Best regards.
Petr