Regresion segmentada:

La regresión segmentada , también conocida como regresión por partes o regresión de varilla rota , es un método en el análisis de regresión en el que la variable independiente se divide en intervalos y un segmento de línea separado se ajusta a cada intervalo. El análisis de regresión segmentado también se puede realizar en datos multivariados dividiendo las diversas variables independientes. La regresión segmentada es útil cuando las variables independientes, agrupadas en diferentes grupos, exhiben diferentes relaciones entre las variables en estas regiones. Los límites entre los segmentos son puntos de interrupción.

Modelo de regresion segmentada

Modelo encontrado en Rpubs: https://rpubs.com/MarkusLoew/12164

Los datos son tomados de https://cran.r-project.org/web/packages/segmented/segmented.pdf

set.seed(10)
x<-1:100
z<-sort(runif(100)) 
y<-sort(2+1.5*pmax(x-35,0)-1.5*pmax(x-70,0)+10*pmax(z-.5,0)+rnorm(100,0,2))

df<-data.frame(z,y)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.5
datos<- ggplot(df, aes(x = y, y =z)) + geom_line() + geom_vline (xintercept = 6) + geom_vline(xintercept = 52.5); datos

Modelo lineal simple

mls = lm(z ~ y, data = df)
summary(mls)
## 
## Call:
## lm(formula = z ~ y, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.14831 -0.07641  0.01514  0.06549  0.17602 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.1586211  0.0127151   12.47   <2e-16 ***
## y           0.0101935  0.0003438   29.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08256 on 98 degrees of freedom
## Multiple R-squared:  0.8997, Adjusted R-squared:  0.8987 
## F-statistic:   879 on 1 and 98 DF,  p-value: < 2.2e-16
coeficiente = coef(mls)
datos + geom_abline(intercept = coeficiente[1],
                    slope = coeficiente[2])

Modelo lineal segmentado : Analisis de puntos de quiebre

library(segmented)
## Warning: package 'segmented' was built under R version 4.0.5
mlseg<- segmented(mls, 
                  seg.Z = ~ y,
                  psi = list(y = c(6,52.5)))
summary(mlseg)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = mls, seg.Z = ~y, psi = list(y = c(6, 52.5)))
## 
## Estimated Break-Point(s):
##           Est. St.Err
## psi1.y  6.012  0.168
## psi2.y 52.869  0.168
## 
## Meaningful coefficients of the linear terms:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.063888   0.003607   17.71   <2e-16 ***
## y            0.048319   0.001289   37.47   <2e-16 ***
## U1.y        -0.043927   0.001303  -33.71       NA    
## U2.y         0.042862   0.001053   40.70       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01431 on 94 degrees of freedom
## Multiple R-Squared: 0.9971,  Adjusted R-squared: 0.997 
## 
## Convergence attained in 2 iter. (rel. change 0)
# Slope 
slope(mlseg)
## $y
##             Est.    St.Err. t value CI(95%).l CI(95%).u
## slope1 0.0483190 0.00128940  37.474 0.0457590 0.0508800
## slope2 0.0043922 0.00018791  23.374 0.0040191 0.0047653
## slope3 0.0472540 0.00103630  45.597 0.0451960 0.0493120
segmentado = fitted(mlseg)
modeloseg = data.frame(y, segmentado)
# Plot 
ggplot(modeloseg, aes(y, segmentado)) + geom_line()

datos + geom_line(data = modeloseg, aes(x = y, y = segmentado), colour = ("green"))

Los datos de clorofila (clf) se miden en el tercio medio, esta es la variable respuesta. El magnesio (Mg) es la variable explicativa para entender el porque la clorofila es mas alta o mas baja.

set.seed(123)

alea = sample(seq(600), 600, F)

cl_enf = rnorm(120, 350, 25) #Clorofila en plantas enfermas
cl_san = rnorm(480, 500, 30) #Clorofila en plantas sanas


estado = rep(c('enferma', 'sana'), c(120, 480))[alea]
cl = c(cl_enf, cl_san)[alea]
mg = sort(runif(600, 150, 250))[alea] #Valores de magnesio 



coor = expand.grid(x = seq(30), y = seq(20)) #Coordenadas 

df = data.frame(coor, estado, cl, mg)

head(df)
##   x y  estado       cl       mg
## 1 1 1    sana 495.1015 218.0295
## 2 2 1    sana 503.0276 226.3107
## 3 3 1    sana 458.3585 177.2977
## 4 4 1    sana 467.0453 236.9847
## 5 5 1    sana 488.1261 180.9994
## 6 6 1 enferma 364.6684 168.4360

Datos

set.seed(123)
Cl = sort(rnorm(600, 400, 20))
mg = sort(runif(600, 150, 250))
estado = as.factor(ifelse(Cl>400, 'sana', 'enferma'))

plot(mg, Cl, col = estado)

+ A medida que sube el Mg tambien hay mayor cantidad de Clorofila

Modelo regresion simple

mod = lm(Cl ~ mg)
summary(mod)
## 
## Call:
## lm(formula = Cl ~ mg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.152  -2.382  -0.045   2.472  32.375 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.710e+02  1.187e+00   228.2   <2e-16 ***
## mg          6.459e-01  5.863e-03   110.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.194 on 598 degrees of freedom
## Multiple R-squared:  0.953,  Adjusted R-squared:  0.953 
## F-statistic: 1.214e+04 on 1 and 598 DF,  p-value: < 2.2e-16
plot(mg, Cl,col = estado, pch = 16, cex = 0.8)
lines(mg, mod$fitted.values, col = 'green', lwd = 2)

Modelo regresion segmentada : Usando la libreria ´segmented´

library(segmented)


mod.seg <- segmented(mod, 
                    seg.Z = ~ mg, psi = c(200))

summary(mod.seg)
## 
##  ***Regression Model with Segmented Relationship(s)***
## 
## Call: 
## segmented.lm(obj = mod, seg.Z = ~mg, psi = c(200))
## 
## 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)
plot(mg, Cl,col = estado, pch = 16, cex = 0.8)
lines(mg, mod.seg$fitted.values, col = 'purple', lwd = 2)

I. Moran

#Extrayendo los residuales 
res = mod.seg$residuals

#Matriz de pesos 
library(ape)
## Warning: package 'ape' was built under R version 4.0.5
d = as.matrix(dist(coor))
di = 1/d
diag(di) = 0

W = di


#Indice de moran
Moran.I(res, W)
## $observed
## [1] 0.06861005
## 
## $expected
## [1] -0.001669449
## 
## $sd
## [1] 0.002358153
## 
## $p.value
## [1] 0

el p-valor del indice de Moran indica que los datos presentan dependencia espacial debido a que no se han involucrado tecnicas espaciales que sirvan para eliminar la dependencia que estan presentando los datos