set.seed(123)
clf = sort(rnorm(600, 400, 20))
Mg = sort(runif(600, 150, 250))
estado = as.factor(ifelse(clf>400, 'sana', 'enferma'))
#clf <- as.factor(clf)
#Mg <- as.factor(Mg)
datos <- data.frame(clf,Mg)
plot(Mg, clf, col = estado)
library(segmented)
model_rseg <- lm(clf ~ Mg, data = datos)
piece <- segmented(model_rseg, seg.Z = ~Mg)
summary(piece)
##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = model_rseg, seg.Z = ~Mg)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.Mg 242.46 0.387
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.791e+02 1.003e+00 278.36 <2e-16 ***
## Mg 6.016e-01 5.068e-03 118.71 <2e-16 ***
## U1.Mg 2.368e+00 1.944e-01 12.18 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.169 on 596 degrees of freedom
## Multiple R-Squared: 0.9733, Adjusted R-squared: 0.9731
##
## Convergence attained in 2 iter. (rel. change 0)
slope(piece)
## $Mg
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.60158 0.0050675 118.710 0.59162 0.61153
## slope2 2.96940 0.1942900 15.283 2.58780 3.35090
intercept(piece)
## $Mg
## Est.
## intercept1 279.08
## intercept2 -295.01
$ y_1= 0.60158x + 279.08$ $y_2 = 2.9694x -295.01 $
my.fit <- fitted(piece)
mymodelo <- data.frame(Mg = datos$Mg, clf = my.fit)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
ggplot(mymodelo, aes(x=Mg, y=clf)) + geom_point(alpha = 0.05) +
geom_point(data = datos, aes(x=Mg, y=clf)) +
geom_abline(intercept = 279.08, slope = 0.60158 , col = "blue") +
geom_abline(intercept = -295.01, slope = 2.96940 , col = "blue") + ylim(c(340,465))