Datos

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)

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


estado = rep(c('enferma', 'sana'), c(120, 480))[alea]
clf = c(clf_enf, clf_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, clf, Mg)

head(df)
##   x y  estado      clf       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

Mapeo de los valores de clorofila

library(ggplot2)

ggplot(data = df)+
  aes(x,y,fill = clf)+
  geom_tile()

Creando los datos ordenados

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

plot(Mg, Clf, col = estado)

Modelo de regresion lineal simple

\[y = \beta_0 + \beta_1x\]

\(\beta_0\) es el intercepto, \(\beta_1\) es la pendiente y
\(x\) es la variable explicativa que en este caso es el Mg

mod = lm(Clf ~ Mg)
summary(mod)
## 
## Call:
## lm(formula = Clf ~ 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, Clf,col = estado, pch = 16, cex = 0.8)
lines(Mg, mod$fitted.values, col = 'blue', lwd = 2)

Modelo de regresion segmentada

La libreria ‘segmented’ permite estimar modelos en casos donde se tienen una o más relaciones segmentadas. Con la funcion ‘segmented’ se pueden estimar un nuevo modelo de regresión que tiene relaciones segmentadas entre la respuesta y una o más variables explicativas y se proporcionan estimaciones de puntos de corte de manera simultanea. El número de puntos de corte de cada relación segmentada se fija mediante el argumento psi, donde se deben especificar los valores iniciales para los puntos de corte, en este caso y la variable segmentada en el argumento seg.Z.

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, Clf,col = estado, pch = 16, cex = 0.8)
lines(Mg, mod.seg$fitted.values, col = 'red', lwd = 2)

Indice de Moran

#Extrayendo los residuales 
res = mod.seg$residuals

#Matriz de pesos 
library(ape)

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