The R Book
Sergio Rocha & Kevin Quiroga
17/11/2020
En el siguiente trabajo se hace una síntesis de los capítulos 14, 16 y 26 del Libro R (\(The~R~Book\), second edition) por \(M.~Crawley\)
Para efectos del aprendizaje los ejemplos mostrados en el libro se han modificado con un enfoque agronómico.
Cap 14 - Count Data - Recuento de Datos
Con el recuento de datos el número \(0\) aparece frecuentemente como valor de la variable respuesta. En el capítulo se tratan los datos como frecuencia (cuantos veces sucede algo), pero no se sabe cuantas veces no sucede. Por ende se utlizan proporciones, donde se sabe el número de veces que sucede algo en específico como las vaces que no sucede.
14.1 - Una regresión con errores de Poisson
Se presenta en el siguiente ejemplo como recuento (número de parcelas con plantas de cebolla podridas por \(Enterobacter~cloacae\) como variable respuesta y la distancia (en kilometros) del rio Bogotá como variable explicativa continua única.
http://www.scielo.org.co/pdf/nova/v16n29/1794-2470-nova-16-29-00071.pdf
clusters <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/clusters.txt")
attach(clusters)
names(clusters)## [1] "Decay" "Distance"Se puede observar una tendencia inversamente proporcional, entre mayor sea la distancia menor seran la pudrición en las plantas de cebolla. Pero, ¿Será esta tendencia significativa? Se debe hacer una regresion de casos vs. distancia con la función glm (con error de Poisson).
model1 <- glm(Decay~Distance,poisson)
summary(model1)##
## Call:
## glm(formula = Decay ~ Distance, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5504 -1.3491 -1.1553 0.3877 3.1304
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.186865 0.188728 0.990 0.3221
## Distance -0.006138 0.003667 -1.674 0.0941 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 149.48 on 93 degrees of freedom
## Residual deviance: 146.64 on 92 degrees of freedom
## AIC: 262.41
##
## Number of Fisher Scoring iterations: 5Si se observa el P valor, se puede evidenciar que la tendencia no es significativa, pero se debe analizar la desviación residual. Bajo el error de Poisson, se asume que la desviación residual debe ser igual a los grados de libertad de los residuales (ya que la varianza y la media deben ser iguales). El hecho de que la desviación residual sea más grande que los grados de libertad de los residuales, indica que se tiene una sobredisperción (extra, inexplicable variación en la respuesta). Se compensa la sobredisperción ajustando el modelo a quasi-Poisson.
model2 <- glm(Decay~Distance,quasipoisson)
summary(model2)##
## Call:
## glm(formula = Decay ~ Distance, family = quasipoisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5504 -1.3491 -1.1553 0.3877 3.1304
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.186865 0.235364 0.794 0.429
## Distance -0.006138 0.004573 -1.342 0.183
##
## (Dispersion parameter for quasipoisson family taken to be 1.555271)
##
## Null deviance: 149.48 on 93 degrees of freedom
## Residual deviance: 146.64 on 92 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Como se observa, al compensar la sobredisperción, esta incrementa el \(p\) valor a 0.183, entonces no hay evidencia convincente que soporte la existencia de una tendencia en la incidencia de la pudrición en las plantas de cebolla con la distancia del río Bogotá. Para ilustrar el modelo ajustado a través de la información, se necesita entender que la función glm con el error Poisson usa el link logaritmico, entonces el parametro estimado y las predicciones desde el modelo (el “predictor lineal”) esta en logaritmos, y necesita ser antilogartimico \(exp(yv)\) antes de dibujar la línea ajustada.
xv <- seq(0,100)
yv <- predict(model2,list(Distance=xv))
y <- c(xv,exp(yv))
fig_dist <- plot_ly(data = clusters, x = ~Distance, y = ~Decay,
marker = list(size = 5,
color = 'brown',
line = list(color = 'brown1',
width = 2)))
fig_dist <- fig_dist %>% layout(title = 'Pudrición vs Distancia',
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE))
#♣fig_dist <- fig_dist %>% add_trace(y = ~y,mode='lines')
fig_dist## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.#lines(xv,exp(yv),col="blue")Note lo extraño que se ve el diagrama de dispersión cuando se tienen los datos de conteo como respuesta. Los valores de la variable respuesta estan en las columnas (porque son todos numeros enteros).
14.2 - Análisis de desviación con recuento de datos
En el siguiente ejemplo, la respuesta variable es un recuento de hojas infectadas de papa con \(Alternaria~solani\), tomadas al azar de un cultivo en la sabana Cundinamarques. Las variables explicativas son:
Patogeno (logico: verdadero o falso)
Etapa fenológica (Tres niveles: joven, medio, viejo)
Tratamiento presente (Si o no)
Rendimiento (Bajo, no normal y normal)
count <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/leaves.txt")
attach(count)
names(count)## [1] "leaves" "Pathogen" "Stage" "Var" "Yield"Siempre es una buena idea con los datos de recuento tener una idea de la distribución de frecuencia general de los recuentos utilizando la función Table:
table(leaves)## leaves
## 0 1 2 3 4 5 6 7
## 314 75 50 32 18 13 7 2La mayoría de plantas (314) no muestran daño en plantas, y el máximo de 7 fue observado en solo dos plantas.
Se empieza la inspección de la información por tabulación el principal efecto de las medias:
tapply(leaves,Pathogen,mean)## FALSE TRUE
## 0.5478723 1.9111111tapply(leaves,Yield, mean)## low normal notnormal
## 1.2814371 0.5833333 0.9357143tapply(leaves,Var,mean)## red white
## 0.6584507 1.2202643tapply(leaves,Stage,mean)## mid old young
## 0.8676471 0.7835821 1.2710280Se ve como si las plantas con patógeno substancialmente tiene un recuento de media más grande que los que no lo tienen; Parcelas con rendimiento no normal y bajo tienen un recuento de media más alto que el normal;Se necesita probar si alguna de estas diferencias son significativas y para evaluar si existen interacciones entre las variables explicativas.
model3<-glm(leaves~Pathogen*Var*Stage*Yield,poisson)
summary(model3)##
## Call:
## glm(formula = leaves ~ Pathogen * Var * Stage * Yield, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.236 -1.022 -0.851 0.520 3.760
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error z value
## (Intercept) 0.4055 0.5774 0.702
## PathogenTRUE 5.8062 2.0520 2.830
## Varwhite -0.9760 0.6405 -1.524
## Stageold -1.1110 0.5986 -1.856
## Stageyoung -1.0986 1.1547 -0.951
## Yieldnormal -1.2384 0.6583 -1.881
## Yieldnotnormal -1.7918 1.1547 -1.552
## PathogenTRUE:Varwhite -3.8247 2.0272 -1.887
## PathogenTRUE:Stageold -3.3089 2.0982 -1.577
## PathogenTRUE:Stageyoung -3.5036 1.7732 -1.976
## Varwhite:Stageold 1.9560 0.6894 2.837
## Varwhite:Stageyoung 2.2569 1.2334 1.830
## PathogenTRUE:Yieldnormal -5.9849 2.1356 -2.802
## PathogenTRUE:Yieldnotnormal -3.7268 1.6464 -2.264
## Varwhite:Yieldnormal 1.1583 0.7704 1.503
## Varwhite:Yieldnotnormal 1.9568 1.2172 1.608
## Stageold:Yieldnormal 0.9280 0.7113 1.305
## Stageyoung:Yieldnormal 1.2384 1.2555 0.986
## Stageold:Yieldnotnormal 1.9293 1.1780 1.638
## Stageyoung:Yieldnotnormal 1.7918 1.6833 1.064
## PathogenTRUE:Varwhite:Stageold 2.6334 2.0906 1.260
## PathogenTRUE:Varwhite:Stageyoung 1.5911 1.7014 0.935
## PathogenTRUE:Varwhite:Yieldnormal 4.6540 2.1470 2.168
## PathogenTRUE:Varwhite:Yieldnotnormal 2.1507 1.4848 1.448
## PathogenTRUE:Stageold:Yieldnormal 0.7871 0.8779 0.897
## PathogenTRUE:Stageyoung:Yieldnormal 4.6268 1.9360 2.390
## PathogenTRUE:Stageold:Yieldnotnormal 2.7137 1.7412 1.559
## PathogenTRUE:Stageyoung:Yieldnotnormal 3.0112 1.0482 2.873
## Varwhite:Stageold:Yieldnormal -2.2210 0.9779 -2.271
## Varwhite:Stageyoung:Yieldnormal -2.5346 1.4309 -1.771
## Varwhite:Stageold:Yieldnotnormal -3.2852 1.4409 -2.280
## Varwhite:Stageyoung:Yieldnotnormal -3.6432 1.8510 -1.968
## PathogenTRUE:Varwhite:Stageold:Yieldnormal NA NA NA
## PathogenTRUE:Varwhite:Stageyoung:Yieldnormal -2.4160 1.9712 -1.226
## PathogenTRUE:Varwhite:Stageold:Yieldnotnormal -0.7991 1.7691 -0.452
## PathogenTRUE:Varwhite:Stageyoung:Yieldnotnormal NA NA NA
## Pr(>|z|)
## (Intercept) 0.48250
## PathogenTRUE 0.00466 **
## Varwhite 0.12756
## Stageold 0.06345 .
## Stageyoung 0.34139
## Yieldnormal 0.05994 .
## Yieldnotnormal 0.12073
## PathogenTRUE:Varwhite 0.05920 .
## PathogenTRUE:Stageold 0.11479
## PathogenTRUE:Stageyoung 0.04816 *
## Varwhite:Stageold 0.00455 **
## Varwhite:Stageyoung 0.06728 .
## PathogenTRUE:Yieldnormal 0.00507 **
## PathogenTRUE:Yieldnotnormal 0.02360 *
## Varwhite:Yieldnormal 0.13272
## Varwhite:Yieldnotnormal 0.10792
## Stageold:Yieldnormal 0.19200
## Stageyoung:Yieldnormal 0.32394
## Stageold:Yieldnotnormal 0.10147
## Stageyoung:Yieldnotnormal 0.28712
## PathogenTRUE:Varwhite:Stageold 0.20780
## PathogenTRUE:Varwhite:Stageyoung 0.34970
## PathogenTRUE:Varwhite:Yieldnormal 0.03018 *
## PathogenTRUE:Varwhite:Yieldnotnormal 0.14749
## PathogenTRUE:Stageold:Yieldnormal 0.36997
## PathogenTRUE:Stageyoung:Yieldnormal 0.01685 *
## PathogenTRUE:Stageold:Yieldnotnormal 0.11911
## PathogenTRUE:Stageyoung:Yieldnotnormal 0.00407 **
## Varwhite:Stageold:Yieldnormal 0.02313 *
## Varwhite:Stageyoung:Yieldnormal 0.07651 .
## Varwhite:Stageold:Yieldnotnormal 0.02261 *
## Varwhite:Stageyoung:Yieldnotnormal 0.04904 *
## PathogenTRUE:Varwhite:Stageold:Yieldnormal NA
## PathogenTRUE:Varwhite:Stageyoung:Yieldnormal 0.22033
## PathogenTRUE:Varwhite:Stageold:Yieldnotnormal 0.65147
## PathogenTRUE:Varwhite:Stageyoung:Yieldnotnormal NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1052.95 on 510 degrees of freedom
## Residual deviance: 736.33 on 477 degrees of freedom
## AIC: 1318
##
## Number of Fisher Scoring iterations: 6Como la desviación residua (736.33) es mucho más grande que los grados de libertad (477), indicando sobredisperción; antes de hacer alguna interpretación se debe ajustar el modelo usando el error quasi-Poisson.
model4<-glm(leaves~Pathogen*Var*Stage*Yield,quasipoisson)
summary(model4)##
## Call:
## glm(formula = leaves ~ Pathogen * Var * Stage * Yield, family = quasipoisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.236 -1.022 -0.851 0.520 3.760
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value
## (Intercept) 0.4055 0.7863 0.516
## PathogenTRUE 5.8062 2.7947 2.078
## Varwhite -0.9760 0.8723 -1.119
## Stageold -1.1110 0.8153 -1.363
## Stageyoung -1.0986 1.5726 -0.699
## Yieldnormal -1.2384 0.8965 -1.381
## Yieldnotnormal -1.7918 1.5726 -1.139
## PathogenTRUE:Varwhite -3.8247 2.7609 -1.385
## PathogenTRUE:Stageold -3.3089 2.8576 -1.158
## PathogenTRUE:Stageyoung -3.5036 2.4149 -1.451
## Varwhite:Stageold 1.9560 0.9389 2.083
## Varwhite:Stageyoung 2.2569 1.6798 1.344
## PathogenTRUE:Yieldnormal -5.9849 2.9085 -2.058
## PathogenTRUE:Yieldnotnormal -3.7268 2.2423 -1.662
## Varwhite:Yieldnormal 1.1583 1.0493 1.104
## Varwhite:Yieldnotnormal 1.9568 1.6578 1.180
## Stageold:Yieldnormal 0.9280 0.9687 0.958
## Stageyoung:Yieldnormal 1.2384 1.7098 0.724
## Stageold:Yieldnotnormal 1.9293 1.6044 1.203
## Stageyoung:Yieldnotnormal 1.7918 2.2924 0.782
## PathogenTRUE:Varwhite:Stageold 2.6334 2.8472 0.925
## PathogenTRUE:Varwhite:Stageyoung 1.5911 2.3171 0.687
## PathogenTRUE:Varwhite:Yieldnormal 4.6540 2.9240 1.592
## PathogenTRUE:Varwhite:Yieldnotnormal 2.1507 2.0222 1.064
## PathogenTRUE:Stageold:Yieldnormal 0.7871 1.1956 0.658
## PathogenTRUE:Stageyoung:Yieldnormal 4.6268 2.6367 1.755
## PathogenTRUE:Stageold:Yieldnotnormal 2.7137 2.3714 1.144
## PathogenTRUE:Stageyoung:Yieldnotnormal 3.0112 1.4275 2.109
## Varwhite:Stageold:Yieldnormal -2.2210 1.3318 -1.668
## Varwhite:Stageyoung:Yieldnormal -2.5346 1.9488 -1.301
## Varwhite:Stageold:Yieldnotnormal -3.2852 1.9623 -1.674
## Varwhite:Stageyoung:Yieldnotnormal -3.6432 2.5209 -1.445
## PathogenTRUE:Varwhite:Stageold:Yieldnormal NA NA NA
## PathogenTRUE:Varwhite:Stageyoung:Yieldnormal -2.4160 2.6846 -0.900
## PathogenTRUE:Varwhite:Stageold:Yieldnotnormal -0.7991 2.4094 -0.332
## PathogenTRUE:Varwhite:Stageyoung:Yieldnotnormal NA NA NA
## Pr(>|t|)
## (Intercept) 0.6063
## PathogenTRUE 0.0383 *
## Varwhite 0.2638
## Stageold 0.1736
## Stageyoung 0.4851
## Yieldnormal 0.1678
## Yieldnotnormal 0.2551
## PathogenTRUE:Varwhite 0.1666
## PathogenTRUE:Stageold 0.2475
## PathogenTRUE:Stageyoung 0.1475
## Varwhite:Stageold 0.0378 *
## Varwhite:Stageyoung 0.1797
## PathogenTRUE:Yieldnormal 0.0402 *
## PathogenTRUE:Yieldnotnormal 0.0972 .
## Varwhite:Yieldnormal 0.2702
## Varwhite:Yieldnotnormal 0.2384
## Stageold:Yieldnormal 0.3386
## Stageyoung:Yieldnormal 0.4693
## Stageold:Yieldnotnormal 0.2297
## Stageyoung:Yieldnotnormal 0.4348
## PathogenTRUE:Varwhite:Stageold 0.3555
## PathogenTRUE:Varwhite:Stageyoung 0.4926
## PathogenTRUE:Varwhite:Yieldnormal 0.1121
## PathogenTRUE:Varwhite:Yieldnotnormal 0.2881
## PathogenTRUE:Stageold:Yieldnormal 0.5107
## PathogenTRUE:Stageyoung:Yieldnormal 0.0799 .
## PathogenTRUE:Stageold:Yieldnotnormal 0.2531
## PathogenTRUE:Stageyoung:Yieldnotnormal 0.0354 *
## Varwhite:Stageold:Yieldnormal 0.0960 .
## Varwhite:Stageyoung:Yieldnormal 0.1940
## Varwhite:Stageold:Yieldnotnormal 0.0948 .
## Varwhite:Stageyoung:Yieldnotnormal 0.1491
## PathogenTRUE:Varwhite:Stageold:Yieldnormal NA
## PathogenTRUE:Varwhite:Stageyoung:Yieldnormal 0.3686
## PathogenTRUE:Varwhite:Stageold:Yieldnotnormal 0.7403
## PathogenTRUE:Varwhite:Stageyoung:Yieldnotnormal NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.854815)
##
## Null deviance: 1052.95 on 510 degrees of freedom
## Residual deviance: 736.33 on 477 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6Lo que primero se debe entender qué son los \(NAs\) en la tabla de medias. No hay información en la tabla de datos con lo cual se pueda estimar esta interacción. Hay aparentemente una interacciónde tres vias entre patogeno, etapa y bajo rendimeinto (\(p~value~=~0.0424\)). Como hay pocos datos para hacer la interacción de 4 vias, se hace una simplificación del modelo, removiendo la interacción de orden más alto.
model5 <- update(model4, ~. -Pathogen:Var:Stage:Yield)
summary(model5)##
## Call:
## glm(formula = leaves ~ Pathogen + Var + Stage + Yield + Pathogen:Var +
## Pathogen:Stage + Var:Stage + Pathogen:Yield + Var:Yield +
## Stage:Yield + Pathogen:Var:Stage + Pathogen:Var:Yield + Pathogen:Stage:Yield +
## Var:Stage:Yield, family = quasipoisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2442 -1.0477 -0.8921 0.5195 3.7613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.4055 0.7849 0.517 0.60567
## PathogenTRUE 3.9549 1.5299 2.585 0.01003 *
## Varwhite -0.9760 0.8707 -1.121 0.26288
## Stageold -1.1121 0.8135 -1.367 0.17228
## Stageyoung -1.0577 1.3450 -0.786 0.43206
## Yieldnormal -1.3027 0.8983 -1.450 0.14768
## Yieldnotnormal -1.3077 1.2453 -1.050 0.29420
## PathogenTRUE:Varwhite -1.9733 1.4674 -1.345 0.17935
## PathogenTRUE:Stageold -1.4495 1.5215 -0.953 0.34121
## PathogenTRUE:Stageyoung -1.6960 1.3007 -1.304 0.19287
## Varwhite:Stageold 1.9587 0.9357 2.093 0.03684 *
## Varwhite:Stageyoung 2.2113 1.4194 1.558 0.11991
## PathogenTRUE:Yieldnormal -3.9246 1.4755 -2.660 0.00808 **
## PathogenTRUE:Yieldnotnormal -2.7324 1.3488 -2.026 0.04334 *
## Varwhite:Yieldnormal 1.2732 1.0407 1.223 0.22178
## Varwhite:Yieldnotnormal 1.4273 1.3336 1.070 0.28505
## Stageold:Yieldnormal 0.9934 0.9705 1.024 0.30656
## Stageyoung:Yieldnormal 1.3469 1.4595 0.923 0.35653
## Stageold:Yieldnotnormal 1.4476 1.2793 1.132 0.25838
## Stageyoung:Yieldnotnormal 0.8630 1.7242 0.500 0.61696
## PathogenTRUE:Varwhite:Stageold 0.7711 1.4515 0.531 0.59549
## PathogenTRUE:Varwhite:Stageyoung -0.2103 1.1404 -0.184 0.85376
## PathogenTRUE:Varwhite:Yieldnormal 2.5007 1.3699 1.825 0.06857 .
## PathogenTRUE:Varwhite:Yieldnotnormal 1.3904 1.0848 1.282 0.20055
## PathogenTRUE:Stageold:Yieldnormal 0.8830 1.1879 0.743 0.45766
## PathogenTRUE:Stageyoung:Yieldnormal 2.4533 1.0471 2.343 0.01954 *
## PathogenTRUE:Stageold:Yieldnotnormal 1.7060 1.2584 1.356 0.17585
## PathogenTRUE:Stageyoung:Yieldnotnormal 2.4941 1.2288 2.030 0.04293 *
## Varwhite:Stageold:Yieldnormal -2.3386 1.3248 -1.765 0.07816 .
## Varwhite:Stageyoung:Yieldnormal -2.8220 1.6238 -1.738 0.08288 .
## Varwhite:Stageold:Yieldnotnormal -2.7807 1.5825 -1.757 0.07953 .
## Varwhite:Stageyoung:Yieldnotnormal -2.4642 1.7181 -1.434 0.15215
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.847991)
##
## Null deviance: 1052.95 on 510 degrees of freedom
## Residual deviance: 737.87 on 479 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6newYd <- Yield
newYd <- as.factor(newYd)
levels(newYd)[c(1,3)] <- "not"
levels(newYd)## [1] "not" "normal"model15 <- glm(leaves~Pathogen+newYd+Pathogen:newYd,
family = quasipoisson)
summary(model15)##
## Call:
## glm(formula = leaves ~ Pathogen + newYd + Pathogen:newYd, family = quasipoisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.3452 -1.1185 -0.9148 0.4352 3.4694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4692 0.1124 -4.173 3.54e-05 ***
## PathogenTRUE 1.4808 0.1484 9.977 < 2e-16 ***
## newYdnormal -0.4021 0.2101 -1.914 0.0562 .
## PathogenTRUE:newYdnormal -0.6583 0.2907 -2.265 0.0239 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 1.858175)
##
## Null deviance: 1052.95 on 510 degrees of freedom
## Residual deviance: 809.41 on 507 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 6Este modelo muestra una interacción muy significativa entre el patógeno y el rendimiento en la determinación del número de hojas dañadas, pero no hay efectos convincentes de la estapa o la variedad. En un caso como este, es útil producir una tabla resumen para resaltar los efectos:
mod_1 <- tapply(leaves,list(Pathogen,Yield),mean)
mod_1## low normal notnormal
## FALSE 0.6893939 0.4184397 0.5436893
## TRUE 3.5142857 0.9523810 2.0270270Este modelo muestra la alta interacción significativa entre el patógeno y el rendimiento; el patógeno agrega una media de alrededor de 0.5 hojas dañadas para parecas con rendimiento normal, pero agrega 2.8 hojas dañadas para parcelas con rendimiento bajo.
Es sencillo convertir la tabla de resumen en un diagrama de barras:
ggplot(data=df, aes(x=Yield, y=Damage_Leaves, fill=Pathogen)) +
geom_bar(stat="identity", position=position_dodge())+
scale_fill_brewer(palette="Greens")+
theme_minimal()
14.3 - Análisis de covarianza con recuento de datos
En el siguiente ejemplo la respuesta es un recuento de la cantidad de especies de plantas en parcelas que tienen diferente biomasa (una variable explicativa continua) y la diferencia del pH del suelo (una vaiable categórica con 3 niveles: alto, medio y bajo)
species <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/species.txt")
attach(species)
names(species)## [1] "pH" "Biomass" "Species"fig_species <- plot_ly(data = species, x = ~Biomass, y = ~Species, color = ~pH)
fig_species## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-modemodel7<-glm(Species~ Biomass*pH,poisson)
summary(model7)##
## Call:
## glm(formula = Species ~ Biomass * pH, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4978 -0.7485 -0.0402 0.5575 3.2297
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.76812 0.06153 61.240 < 2e-16 ***
## Biomass -0.10713 0.01249 -8.577 < 2e-16 ***
## pHlow -0.81557 0.10284 -7.931 2.18e-15 ***
## pHmid -0.33146 0.09217 -3.596 0.000323 ***
## Biomass:pHlow -0.15503 0.04003 -3.873 0.000108 ***
## Biomass:pHmid -0.03189 0.02308 -1.382 0.166954
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 452.346 on 89 degrees of freedom
## Residual deviance: 83.201 on 84 degrees of freedom
## AIC: 514.39
##
## Number of Fisher Scoring iterations: 4Como se evidencia no se tiene sobredispersión (grados de libertad \(>\) desviación residual). Se puede probar la necesidad de diferentes pendientes comparando este modelo máximo (con seis parámetros) con un modelo más simple con diferentes intersecciones pero la misma pendiente (cuatro parámetros):
model8 <- glm(Species~Biomass+pH,poisson)
summary(model8)##
## Call:
## glm(formula = Species ~ Biomass + pH, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5959 -0.6989 -0.0737 0.6647 3.5604
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.84894 0.05281 72.885 < 2e-16 ***
## Biomass -0.12756 0.01014 -12.579 < 2e-16 ***
## pHlow -1.13639 0.06720 -16.910 < 2e-16 ***
## pHmid -0.44516 0.05486 -8.114 4.88e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 452.346 on 89 degrees of freedom
## Residual deviance: 99.242 on 86 degrees of freedom
## AIC: 526.43
##
## Number of Fisher Scoring iterations: 4anova(model7,model8,test="Chi")## Analysis of Deviance Table
##
## Model 1: Species ~ Biomass * pH
## Model 2: Species ~ Biomass + pH
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 84 83.201
## 2 86 99.242 -2 -16.04 0.0003288 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Las pendientes son significativamente muy diferentes (p = 0.000 33), por lo está justificado en retener los más complicado del modelo1.
Finalmente, se sibuja las líneas ajustadas a través del diagrama de dispersión, usando \(predict\). Se necesita especificar valores para todas las variables explicativas en el modelo, a saber, biomasa (una variable continua) y pH del suelo (variable categórica con tres niveles). Primero, la variable continua para el eje x:
xv_1 <- seq(0,10,0.1)A continuación, se necesita proporcionar un vector de niveles de factor para el pH del suelo y este vector debe tener exactamente la misma longitud como el vector de valores x (\(longitud~=~101\)). Es más sencillo utilizar los niveles de factor de pH en el orden en que aparecen:
pH <- species$pH
pH <- as.factor(pH)
levels(pH)## [1] "high" "low" "mid"Primero se traza la línea para el pH alto, recordando anular (\(antilog\)) las predicciones:
pHs<-factor(rep("high",101))
yv_1 <- predict(model7,list(Biomass=xv_1,pH=pHs))
pHs_1 <-factor(rep("low",101))
yv_2 <- predict(model7,list(Biomass=xv_1,pH=pHs_1))
pHs_2 <-factor(rep("mid",101))
yv_3 <- predict(model7,list(Biomass=xv_1,pH=pHs_2))
plot(Biomass,Species,type="n")
spp<-split(Species,pH)
bio<-split(Biomass,pH)
points(bio[[1]],spp[[1]],pch=16,col="red")
points(bio[[2]],spp[[2]],pch=16,col="green")
points(bio[[3]],spp[[3]],pch=16,col="blue")
legend(x=8,y=45.5,legend=c("high","low","medium"), pch=c(16,16,16),col=c("red","green","blue"),title="pH")
lines(xv_1,exp(yv_1),col='red')
lines(xv_1,exp(yv_2),col='green')
lines(xv_1,exp(yv_3),col='blue')
14.4 - Distribución de frecuencias
Aquí hay datos sobre el número de quiebras en 80 distritos. La pregunta es si hay alguna evidencia que algunos distritos muestran un número de casos mayor de lo esperado. ¿Qué esperaríamos? Por supuesto que deberíamos esperar alguna variación, pero ¿cuánto exactamente? Bueno, eso depende de nuestro modelo de proceso. Quizás el modelo más simple es que no pasa absolutamente nada, y que cada caso de quiebra individual es absolutamente independientes unos de otros. Eso lleva a la predicción de que el número de casos por distrito seguirá un proceso de Poisson, una distribución en la que la varianza es igual a la media. Observe lo que muestra de datos.
case.book <- read.csv("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/cases.txt", sep="")
attach(case.book)
names(case.book)## [1] "cases"Primero se necesita contar el número de distritos sin casos, un caso, dos casos, etc.
frequencies <- table(cases)
frequencies## cases
## 0 1 2 3 4 5 6 7 8 9 10
## 34 14 10 7 4 5 2 1 1 1 1No hubo ningún caso en 34 distritos, pero un distrito tuvo 10 casos. Una buena forma de proceder es comparar nuestra distribución (llamada frecuencias) con la distribución que se observaría si los datos realmente procedieran de una distribución de Poisson como postula nuestro modelo. Podemos usar la función \(dpois\) para calcular la densidad de probabilidad de cada una de las 11 frecuencias de 0 a 10 (multiplicamos la probabilidad producida por \(dpois\) por la muestra total de 80 para obtener las frecuencias predichas) .Se necesita calcular el número medio de casos por distrito: este es el único parámetro de la distribución de Poisson:
cases <- case.book$cases
cases <- as.vector(cases)
mean(cases)## [1] 1.775El plan es dibujar dos distribuciones una al lado de la otra, por lo que se configura la región de trazado:
fig_poi <- plot_ly(y=~frequencies,type = 'bar')
fig_poi <- fig_poi %>% layout(title='Frecuencias',
xaxis =list(title = "Cases"),
yaxis =list(title = "Frequiencies"))
fig_poifig_poi_1 <- plot_ly(y=~cases_1,type = 'bar')
fig_poi_1 <- fig_poi_1 %>% layout(title='Poisson',
xaxis =list(title = "Cases"),
yaxis =list(title = "Frequiencies"))
fig_poi_1Las distribuciones son muy diferentes: la moda de los datos observados es 0, pero la moda de la distribución de Poisson con la misma media es 1; los datos observados contenían ejemplos de 8, 9 y 10 casos, pero estos serían muy poco probables en un proceso de Poisson. Diríamos que los datos observados están altamente agregados tienen una razón varianza-media mucho mayor que 1 (la distribución de Poisson, por supuesto, tiene una razón varianza-media de 1):
var(cases)/mean(cases)## [1] 2.99483Entonces, si los datos no se distribuyen según Poisson, ¿cómo se distribuyen? Una buena distribución candidata donde la relación varianza-media es así de grande (alrededor de 3,0) es la distribución binomial negativa (consulte la pág. 315). Esta es una distribución de dos parámetros: el primer parámetro es el número medio de casos (1.775) y el segundo se llama parámetro de agrupamiento, k (mide el grado de agregación en los datos: valores pequeños de k (k <1) muestran una alta agregación, mientras que los valores grandes de k (k> 5) muestran aleatoriedad). Podemos obtener una estimación aproximada de la magnitud de k de
\[\hat k = \frac{\hat x}{S^2-\hat x}\]
Se puede trabajar esto afuera
mean(cases)^2/(var(cases)-mean(cases))## [1] 0.8898003Entonces se trabajaá con k = 0.89. ¿Cómo calular las frecuencias esperadas? La función de densidad para la distribución binomial negativa es \(dnbinom\) y tiene tres argumentos: la frecuencia para la que se quiere la probabilidad (en este caso de 0 a 10), el número de éxitos (en este caso 1) y el número medio de casos. (1,775); se multiplica por el número total de casos (80) para obtener las frecuencias esperadas
exp <- dnbinom(0:10,1,mu=1.775)*80Dibujar una sola figura en la que las frecuencias observadas y esperadas se dibujan una al lado de la otra. El truco consiste en producir un nuevo vector (llamado “both”) que sea el doble de largo que los vectores de frecuencia observados y esperados (2 × 11 = 22). Luego, colocar las frecuencias observadas en los elementos impares (usando el módulo 2 para calcular los valores de los subíndices), y las frecuencias esperadas en los elementos pares:
df_both <- data.frame(frequencies,exp)
df_both## cases Freq exp
## 1 0 34 28.8288288
## 2 1 14 18.4400617
## 3 2 10 11.7949944
## 4 3 7 7.5445460
## 5 4 4 4.8257907
## 6 5 5 3.0867670
## 7 6 2 1.9744185
## 8 7 1 1.2629164
## 9 8 1 0.8078114
## 10 9 1 0.5167082
## 11 10 1 0.3305070fig_both <- plot_ly(df_both,y=~frequencies,
type = 'bar',name='Observed')
fig_both <- fig_both %>% add_trace(y=~exp,name='Expected')
fig_both <- fig_both %>% layout(yaxis = list(title = 'Frequency'))
fig_bothEl ajuste a la distribución binomial negativa es mucho mejor que con la distribución de Poisson, especialmente en la cola de la derecha. Pero los datos observados tienen demasiados 0 y muy pocos 1 para ser representados perfectamente por una distribución binomial negativa. Si desea cuantificar la falta de ajuste entre las distribuciones de frecuencia observadas y esperadas, puede calcular Chi cuadro de Pearson \(\sum(O-E)^2/E\) basado en el numero de comparaciones que tiene una frecuencia esperada mayor a 4:
exp## [1] 28.8288288 18.4400617 11.7949944 7.5445460 4.8257907 3.0867670
## [7] 1.9744185 1.2629164 0.8078114 0.5167082 0.3305070Si se acumulan las seis frecuencias más a la derecha, entonces todos los valores de \(exp\) serán mayores que 4. Los grados de libertad vienen dados por el número de comparaciones legítimas (6) menos el número de parámetros estimados a partir de los datos (2 en nuestro caso ) menos 1 (para contingencia, porque la frecuencia total debe sumar 80), es decir, 3 grados de libertad. Reducimos las longitudes de los vectores observados y esperados, creando un intervalo superior llamado 5+ para 5’ o más":
cs <- factor(0:10)
levels(cs)[6:11]<-"5+"
levels(cs)## [1] "0" "1" "2" "3" "4" "5+"Ahora haga dos vestores cortos \(of\) y \(ef\) (para frecuencia \('observada'\) y \('esperada'\))
ef <- as.vector(tapply(exp,cs,sum))
of <- as.vector(tapply(frequencies,cs,sum))Finalmente, se puede calcular el valor de chi-cuadrado midiendo la diferencia entre las distribuciones de frecuencia observadas y esperadas, y usar \(1-pchisq\) para calcular el valor p:
sum((of-ef)^2/ef)## [1] 3.5941451-pchisq(3.594145,3)## [1] 0.3087555Se concluye que una descripción binomial negativa de estos datos es razonable (las distribuciones observadas y esperadas no son significativamente diferentes, p = 0.31).
Sobredispersión en modelos logaritmicos-lineales
Los datos analizados en esta sección se refieren a niños de Guateque – Boyaca - Colombia, que fueron clasificados por sexo (con dos niveles: masculino (M) y femenino (F)), cultura (también con dos niveles: Campesino (A) y no (N)), grupo de edad (con cuatro niveles: F0 (primario), F1, F2 y F3) y estado del alumno (con dos niveles: promedio (AL) y lento (SL)). La variable de respuesta es un recuento del número de días ausentes de la escuela en un año escolar en particular (Días).
data(quine)
attach(quine)
names(quine)## [1] "Eth" "Sex" "Age" "Lrn" "Days"datatable(quine,filter='top',class = 'cell-border stripe', options = list(
pageLength = 20, autoWidth = TRUE))Comienza con un modelo log-lineal para los recuentos y se ajusta un modelo máximo que contiene todos los factores y todas sus interacciones:
model9 <- glm(Days~Eth*Sex*Age*Lrn,poisson)
summary(model9)##
## Call:
## glm(formula = Days ~ Eth * Sex * Age * Lrn, family = poisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.3872 -2.5129 -0.4205 1.7424 6.6783
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0564 0.1085 28.178 < 2e-16 ***
## EthN -0.1386 0.1590 -0.872 0.383394
## SexM -0.4914 0.1648 -2.982 0.002860 **
## AgeF1 -0.6227 0.1712 -3.638 0.000275 ***
## AgeF2 -2.3632 0.7154 -3.303 0.000955 ***
## AgeF3 -0.3784 0.1393 -2.717 0.006592 **
## LrnSL -1.9577 0.5875 -3.333 0.000860 ***
## EthN:SexM -0.7524 0.2682 -2.806 0.005021 **
## EthN:AgeF1 0.1029 0.2408 0.427 0.669209
## EthN:AgeF2 -0.5546 1.2350 -0.449 0.653410
## EthN:AgeF3 0.0633 0.2008 0.315 0.752564
## SexM:AgeF1 0.4092 0.3038 1.347 0.178074
## SexM:AgeF2 3.1098 0.7296 4.262 2.02e-05 ***
## SexM:AgeF3 1.1145 0.2001 5.570 2.55e-08 ***
## EthN:LrnSL 2.2588 0.6314 3.578 0.000347 ***
## SexM:LrnSL 1.5900 0.6305 2.522 0.011673 *
## AgeF1:LrnSL 2.6421 0.6059 4.361 1.30e-05 ***
## AgeF2:LrnSL 4.8585 0.9212 5.274 1.33e-07 ***
## AgeF3:LrnSL NA NA NA NA
## EthN:SexM:AgeF1 -0.3105 0.5432 -0.572 0.567587
## EthN:SexM:AgeF2 0.3469 1.2620 0.275 0.783401
## EthN:SexM:AgeF3 0.8329 0.3122 2.668 0.007627 **
## EthN:SexM:LrnSL -0.1639 0.7024 -0.233 0.815496
## EthN:AgeF1:LrnSL -3.5493 0.6715 -5.286 1.25e-07 ***
## EthN:AgeF2:LrnSL -3.3315 1.3856 -2.404 0.016202 *
## EthN:AgeF3:LrnSL NA NA NA NA
## SexM:AgeF1:LrnSL -2.4285 0.7100 -3.420 0.000626 ***
## SexM:AgeF2:LrnSL -4.1914 0.9555 -4.387 1.15e-05 ***
## SexM:AgeF3:LrnSL NA NA NA NA
## EthN:SexM:AgeF1:LrnSL 2.1711 0.8924 2.433 0.014985 *
## EthN:SexM:AgeF2:LrnSL 2.1029 1.4330 1.467 0.142254
## EthN:SexM:AgeF3:LrnSL NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2073.5 on 145 degrees of freedom
## Residual deviance: 1173.9 on 118 degrees of freedom
## AIC: 1818.4
##
## Number of Fisher Scoring iterations: 5A continuación, verifique la desviación residual para ver si hay sobredispersión. Recuerde que la desviación residual debe ser igual a los grados de libertad residuales si el supuesto de errores de Poisson es apropiado. Aquí es 1173,9 sobre 118 d.f., lo que indica una dispersión excesiva por un factor de aproximadamente 10. Esto es demasiado grande para ignorarlo, por lo que antes de embarcarse en la simplificación del modelo, intente un enfoque diferente, utilizando errores de cuasi-Poisson para explicar la dispersión excesiva:
model10 <- glm(Days~Eth*Sex*Age*Lrn,quasipoisson)
summary(model10)##
## Call:
## glm(formula = Days ~ Eth * Sex * Age * Lrn, family = quasipoisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.3872 -2.5129 -0.4205 1.7424 6.6783
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0564 0.3346 9.135 2.22e-15 ***
## EthN -0.1386 0.4904 -0.283 0.7780
## SexM -0.4914 0.5082 -0.967 0.3356
## AgeF1 -0.6227 0.5281 -1.179 0.2407
## AgeF2 -2.3632 2.2066 -1.071 0.2864
## AgeF3 -0.3784 0.4296 -0.881 0.3802
## LrnSL -1.9577 1.8120 -1.080 0.2822
## EthN:SexM -0.7524 0.8272 -0.910 0.3649
## EthN:AgeF1 0.1029 0.7427 0.139 0.8901
## EthN:AgeF2 -0.5546 3.8094 -0.146 0.8845
## EthN:AgeF3 0.0633 0.6194 0.102 0.9188
## SexM:AgeF1 0.4092 0.9372 0.437 0.6632
## SexM:AgeF2 3.1098 2.2506 1.382 0.1696
## SexM:AgeF3 1.1145 0.6173 1.806 0.0735 .
## EthN:LrnSL 2.2588 1.9474 1.160 0.2484
## SexM:LrnSL 1.5900 1.9448 0.818 0.4152
## AgeF1:LrnSL 2.6421 1.8688 1.414 0.1601
## AgeF2:LrnSL 4.8585 2.8413 1.710 0.0899 .
## AgeF3:LrnSL NA NA NA NA
## EthN:SexM:AgeF1 -0.3105 1.6756 -0.185 0.8533
## EthN:SexM:AgeF2 0.3469 3.8928 0.089 0.9291
## EthN:SexM:AgeF3 0.8329 0.9629 0.865 0.3888
## EthN:SexM:LrnSL -0.1639 2.1666 -0.076 0.9398
## EthN:AgeF1:LrnSL -3.5493 2.0712 -1.714 0.0892 .
## EthN:AgeF2:LrnSL -3.3315 4.2739 -0.779 0.4373
## EthN:AgeF3:LrnSL NA NA NA NA
## SexM:AgeF1:LrnSL -2.4285 2.1901 -1.109 0.2697
## SexM:AgeF2:LrnSL -4.1914 2.9472 -1.422 0.1576
## SexM:AgeF3:LrnSL NA NA NA NA
## EthN:SexM:AgeF1:LrnSL 2.1711 2.7527 0.789 0.4319
## EthN:SexM:AgeF2:LrnSL 2.1029 4.4203 0.476 0.6351
## EthN:SexM:AgeF3:LrnSL NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 9.514226)
##
## Null deviance: 2073.5 on 145 degrees of freedom
## Residual deviance: 1173.9 on 118 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Observe que ciertas interacciones tienen vacios (mostradas por filas con NA). Estos no se han estimado debido a que faltan combinaciones de nivel de factor, como lo indican los ceros en la siguiente tabla:
ftable(table(Eth,Sex,Age,Lrn))## Lrn AL SL
## Eth Sex Age
## A F F0 4 1
## F1 5 10
## F2 1 8
## F3 9 0
## M F0 5 3
## F1 2 3
## F2 7 4
## F3 7 0
## N F F0 4 1
## F1 6 11
## F2 1 9
## F3 10 0
## M F0 6 3
## F1 2 7
## F2 7 3
## F3 7 0La mayor parte de esto ocurre porque los aprendices lentos nunca entran en el Form 3.
Desafortunadamente, el citerio de información de Akaike no está definido para este modelo, por lo que no se puede automatizar la simplificación usando \(step\) o \(stepAIC\). Se necesita hacer la simplificación del modelo a mano, por lo tanto, recordando hacer pruebas F (no chi-cuadrado) debido a la dispersión vertical. Aquí está el último paso de la simplificación antes de obtener el modelo mínimo adecuado. ¿Se necesita la interacción de aprendizaje de la edad?
model11_1 <- glm(Days~Eth*Sex*Age*Lrn,quasipoisson)
model11_2 <- update(model11_1, ~. -Eth:Sex:Age:Lrn)
model11_3 <- update(model11_2, ~. -Eth:Age:Lrn)
model11_4 <- update(model11_3, ~. -Sex:Age:Lrn)
model11 <- update(model11_4, ~. -Eth:Age:Sex)
model12 <- update(model11,~. - Age:Lrn)
anova(model11,model12,test="F")## Analysis of Deviance Table
##
## Model 1: Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age +
## Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn
## Model 2: Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age +
## Eth:Lrn + Sex:Lrn + Eth:Sex:Lrn
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 127 1280.5
## 2 129 1301.1 -2 -20.558 1.0305 0.3598Entonces, aquí está el modelo mínimo adecuado con errores de cuasi-Poisson:
summary(model12)##
## Call:
## glm(formula = Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age +
## Sex:Age + Eth:Lrn + Sex:Lrn + Eth:Sex:Lrn, family = quasipoisson)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -7.1369 -2.6852 -0.5919 1.6040 7.1049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.83161 0.30489 9.287 4.98e-16 ***
## EthN 0.09821 0.38631 0.254 0.79973
## SexM -0.56268 0.38877 -1.447 0.15023
## AgeF1 -0.20878 0.35933 -0.581 0.56224
## AgeF2 0.16223 0.37481 0.433 0.66586
## AgeF3 -0.25584 0.37855 -0.676 0.50036
## LrnSL 0.50311 0.30798 1.634 0.10479
## EthN:SexM -0.24554 0.37347 -0.657 0.51206
## EthN:AgeF1 -0.68742 0.46823 -1.468 0.14450
## EthN:AgeF2 -1.07361 0.42449 -2.529 0.01264 *
## EthN:AgeF3 0.01879 0.42914 0.044 0.96513
## SexM:AgeF1 -0.26358 0.50673 -0.520 0.60385
## SexM:AgeF2 0.94531 0.43530 2.172 0.03171 *
## SexM:AgeF3 1.35285 0.42933 3.151 0.00202 **
## EthN:LrnSL -0.65154 0.45857 -1.421 0.15778
## SexM:LrnSL -0.29570 0.41144 -0.719 0.47363
## EthN:SexM:LrnSL 1.60463 0.57113 2.810 0.00573 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 9.833478)
##
## Null deviance: 2073.5 on 145 degrees of freedom
## Residual deviance: 1301.1 on 129 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5Existe una interacción de tres vías muy significativa entre origen étnico, sexo y dificultad de aprendizaje; Los niños no aborígenes de aprendizaje lento tenían más probabilidades de estar ausentes que los niños no aborígenes sin dificultades de aprendizaje.
ftable(tapply(Days,list(Eth,Sex,Lrn),mean))## AL SL
##
## A F 14.47368 27.36842
## M 22.28571 20.20000
## N F 13.14286 7.00000
## M 13.36364 17.00000Sin embargo, tenga en cuenta que entre los alumnos sin dificultades de aprendizaje son los niños aborígenes los que faltan más días, y son las niñas aborígenes con dificultades de aprendizaje las que tienen la tasa más alta de ausentismo en general.
14.6 - Errores binomiales negativos
En lugar de utilizar errores de cuasi-Poisson (como se indicó anteriormente), se podría utilizar un modelo binomial negativo. Esto está en la biblioteca \(MASS\) e involucra la función \(glm.nb\). El modelado procede exactamente de la misma manera que con un típico GLM:
model.nb1 <- glm.nb(Days~Eth*Sex*Age*Lrn)
summary(model.nb1,cor=F)##
## Call:
## glm.nb(formula = Days ~ Eth * Sex * Age * Lrn, init.theta = 1.928360145,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.2377 -0.9079 -0.2019 0.5173 1.7043
##
## Coefficients: (4 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0564 0.3760 8.128 4.38e-16 ***
## EthN -0.1386 0.5334 -0.260 0.795023
## SexM -0.4914 0.5104 -0.963 0.335653
## AgeF1 -0.6227 0.5125 -1.215 0.224334
## AgeF2 -2.3632 1.0770 -2.194 0.028221 *
## AgeF3 -0.3784 0.4546 -0.832 0.405215
## LrnSL -1.9577 0.9967 -1.964 0.049493 *
## EthN:SexM -0.7524 0.7220 -1.042 0.297400
## EthN:AgeF1 0.1029 0.7123 0.144 0.885175
## EthN:AgeF2 -0.5546 1.6798 -0.330 0.741297
## EthN:AgeF3 0.0633 0.6396 0.099 0.921159
## SexM:AgeF1 0.4092 0.8299 0.493 0.621973
## SexM:AgeF2 3.1098 1.1655 2.668 0.007624 **
## SexM:AgeF3 1.1145 0.6365 1.751 0.079926 .
## EthN:LrnSL 2.2588 1.3019 1.735 0.082743 .
## SexM:LrnSL 1.5900 1.1499 1.383 0.166750
## AgeF1:LrnSL 2.6421 1.0821 2.442 0.014618 *
## AgeF2:LrnSL 4.8585 1.4423 3.369 0.000755 ***
## AgeF3:LrnSL NA NA NA NA
## EthN:SexM:AgeF1 -0.3105 1.2055 -0.258 0.796735
## EthN:SexM:AgeF2 0.3469 1.7965 0.193 0.846875
## EthN:SexM:AgeF3 0.8329 0.8970 0.929 0.353092
## EthN:SexM:LrnSL -0.1639 1.5250 -0.107 0.914411
## EthN:AgeF1:LrnSL -3.5493 1.4270 -2.487 0.012876 *
## EthN:AgeF2:LrnSL -3.3315 2.0919 -1.593 0.111256
## EthN:AgeF3:LrnSL NA NA NA NA
## SexM:AgeF1:LrnSL -2.4285 1.4201 -1.710 0.087246 .
## SexM:AgeF2:LrnSL -4.1914 1.6201 -2.587 0.009679 **
## SexM:AgeF3:LrnSL NA NA NA NA
## EthN:SexM:AgeF1:LrnSL 2.1711 1.9192 1.131 0.257963
## EthN:SexM:AgeF2:LrnSL 2.1029 2.3444 0.897 0.369718
## EthN:SexM:AgeF3:LrnSL NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.9284) family taken to be 1)
##
## Null deviance: 272.29 on 145 degrees of freedom
## Residual deviance: 167.45 on 118 degrees of freedom
## AIC: 1097.3
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.928
## Std. Err.: 0.269
##
## 2 x log-likelihood: -1039.324La salida es ligeramente diferente a la de un GLM convencional: ve el parámetro binomial negativo estimado (aquí llamado Theta, pero conocido como k, e igual a 1.928) y su error estándar aproximado (0.269) y 2 veces la probabilidad logarítmica (contraste esto con la desviación residual de nuestro modelo cuasi-Poisson, que fue 1301.1; ver arriba). Tenga en cuenta que la desviación residual en el modelo binomial negativo (167,45) no es 2 veces la probabilidad logarítmica.
Una ventaja del modelo binomial negativo sobre el cuasi-Poisson es que se puede automatizar la simplificación del modelo con \(stepAIC\):
model.nb2 <- stepAIC(model.nb1)## Start: AIC=1095.32
## Days ~ Eth * Sex * Age * Lrn
##
## Df AIC
## - Eth:Sex:Age:Lrn 2 1092.7
## <none> 1095.3
##
## Step: AIC=1092.73
## Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age +
## Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Age + Eth:Sex:Lrn +
## Eth:Age:Lrn + Sex:Age:Lrn
##
## Df AIC
## - Eth:Sex:Age 3 1089.4
## <none> 1092.7
## - Eth:Sex:Lrn 1 1093.3
## - Eth:Age:Lrn 2 1094.7
## - Sex:Age:Lrn 2 1095.0
##
## Step: AIC=1089.41
## Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age +
## Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn +
## Sex:Age:Lrn
##
## Df AIC
## <none> 1089.4
## - Sex:Age:Lrn 2 1091.1
## - Eth:Age:Lrn 2 1091.2
## - Eth:Sex:Lrn 1 1092.5summary(model.nb2,cor=F)##
## Call:
## glm.nb(formula = Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age +
## Sex:Age + Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn +
## Sex:Age:Lrn, init.theta = 1.865343469, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.1387 -0.9777 -0.2000 0.5349 1.7630
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1693 0.3411 9.292 < 2e-16 ***
## EthN -0.3560 0.4210 -0.845 0.397848
## SexM -0.6920 0.4138 -1.672 0.094459 .
## AgeF1 -0.6405 0.4638 -1.381 0.167329
## AgeF2 -2.4576 0.8675 -2.833 0.004612 **
## AgeF3 -0.5880 0.3973 -1.480 0.138885
## LrnSL -1.0264 0.7378 -1.391 0.164179
## EthN:SexM -0.3562 0.3854 -0.924 0.355364
## EthN:AgeF1 0.1500 0.5644 0.266 0.790400
## EthN:AgeF2 -0.3833 0.5640 -0.680 0.496746
## EthN:AgeF3 0.4719 0.4542 1.039 0.298824
## SexM:AgeF1 0.2985 0.6047 0.494 0.621597
## SexM:AgeF2 3.2904 0.8941 3.680 0.000233 ***
## SexM:AgeF3 1.5412 0.4548 3.389 0.000702 ***
## EthN:LrnSL 0.9651 0.7753 1.245 0.213255
## SexM:LrnSL 0.5457 0.8013 0.681 0.495873
## AgeF1:LrnSL 1.6231 0.8222 1.974 0.048373 *
## AgeF2:LrnSL 3.8321 1.1054 3.467 0.000527 ***
## AgeF3:LrnSL NA NA NA NA
## EthN:SexM:LrnSL 1.3578 0.5914 2.296 0.021684 *
## EthN:AgeF1:LrnSL -2.1013 0.8728 -2.408 0.016058 *
## EthN:AgeF2:LrnSL -1.8260 0.8774 -2.081 0.037426 *
## EthN:AgeF3:LrnSL NA NA NA NA
## SexM:AgeF1:LrnSL -1.1086 0.9409 -1.178 0.238671
## SexM:AgeF2:LrnSL -2.8800 1.1550 -2.493 0.012651 *
## SexM:AgeF3:LrnSL NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8653) family taken to be 1)
##
## Null deviance: 265.27 on 145 degrees of freedom
## Residual deviance: 167.44 on 123 degrees of freedom
## AIC: 1091.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.865
## Std. Err.: 0.258
##
## 2 x log-likelihood: -1043.409Retomar la simplificación del modelo donde lo deja \(AIC\):
model.nb3 <- update(model.nb2,~. - Sex:Age:Lrn)
anova(model.nb3,model.nb2)## Likelihood ratio tests of Negative Binomial Models
##
## Response: Days
## Model
## 1 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age + Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn
## 2 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age + Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn + Sex:Age:Lrn
## theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 1.789507 125 -1049.111
## 2 1.865343 123 -1043.409 1 vs 2 2 5.701942 0.05778817Debido a que no se toma en cuenta, el sexo marginalmente significativo por edad por interacción de aprendizaje no sobrevive a una prueba de eliminación (p = 0.058), ni tampoco por origen étnico por edad por aprendizaje (p = 0.115) ni edad por aprendizaje (p = 0.150)
model.nb4<-update(model.nb3,~. - Eth:Age:Lrn)
anova(model.nb3,model.nb4)## Likelihood ratio tests of Negative Binomial Models
##
## Response: Days
## Model
## 1 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age + Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn
## 2 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age + Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn + Eth:Age:Lrn
## theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 1.724987 127 -1053.431
## 2 1.789507 125 -1049.111 1 vs 2 2 4.320086 0.1153202model.nb5<-update(model.nb4,~. - Age:Lrn)
anova(model.nb4,model.nb5)## Likelihood ratio tests of Negative Binomial Models
##
## Response: Days
## Model
## 1 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age + Eth:Lrn + Sex:Lrn + Eth:Sex:Lrn
## 2 Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age + Sex:Age + Eth:Lrn + Sex:Lrn + Age:Lrn + Eth:Sex:Lrn
## theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 1.678620 129 -1057.219
## 2 1.724987 127 -1053.431 1 vs 2 2 3.787823 0.150482summary(model.nb5,cor=F)##
## Call:
## glm.nb(formula = Days ~ Eth + Sex + Age + Lrn + Eth:Sex + Eth:Age +
## Sex:Age + Eth:Lrn + Sex:Lrn + Eth:Sex:Lrn, init.theta = 1.678619829,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0246 -0.9449 -0.2228 0.4847 1.9002
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.91755 0.32626 8.942 < 2e-16 ***
## EthN 0.05666 0.39515 0.143 0.88598
## SexM -0.55047 0.39014 -1.411 0.15825
## AgeF1 -0.32379 0.38373 -0.844 0.39878
## AgeF2 -0.06383 0.42046 -0.152 0.87933
## AgeF3 -0.34854 0.39128 -0.891 0.37305
## LrnSL 0.57697 0.33382 1.728 0.08392 .
## EthN:SexM -0.41608 0.37491 -1.110 0.26708
## EthN:AgeF1 -0.56613 0.43162 -1.312 0.18965
## EthN:AgeF2 -0.89577 0.42950 -2.086 0.03702 *
## EthN:AgeF3 0.08467 0.44010 0.192 0.84744
## SexM:AgeF1 -0.08459 0.45324 -0.187 0.85195
## SexM:AgeF2 1.13752 0.45192 2.517 0.01183 *
## SexM:AgeF3 1.43124 0.44365 3.226 0.00126 **
## EthN:LrnSL -0.78724 0.43058 -1.828 0.06750 .
## SexM:LrnSL -0.47437 0.45908 -1.033 0.30147
## EthN:SexM:LrnSL 1.75289 0.58341 3.005 0.00266 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.6786) family taken to be 1)
##
## Null deviance: 243.98 on 145 degrees of freedom
## Residual deviance: 168.03 on 129 degrees of freedom
## AIC: 1093.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.679
## Std. Err.: 0.227
##
## 2 x log-likelihood: -1057.219El modelo mínimo adecuado, por lo tanto, contiene exactamente los mismos términos que se obtuvieron con cuasi-Poisson, pero los niveles de significancia son más altos (por ejemplo, la interacción de tres vías tiene p = 0.002 66 en comparación con p = 0.005 73). es necesario trazar el modelo para verificar los supuestos:
plot(model.nb5)
La varianza se comporta bien y los residuos están distribuidos casi normalmente. La combinación de valores p bajos y la capacidad de utilizar stepAIC hace que glm.nb sea una función de modelado muy útil para datos de recuento como estos.
Cap 16 - Proporción de información
Una clase importante de problemas implica contar datos en proporciones como:
Estudios sobre tasas de mortalidad.
Tasas de infección de enfermedades.
Respuestas a cuestionarios.
Proporción que responde al tratamiento clínico.
Proporción que admite intenciones de voto particulares.
Proporciones de sexos.
Datos sobre la respuesta proporcional a un tratamiento experimental.
Lo que todos estos tienen en común es que se sabe cuántos de los objetos experimentales están en una categoría (muertos, insolventes, masculinos o infectados) y también se sabe cuántos hay en otra (vivos, solventes, femeninos o no infectados). Esto contrasta con los datos del recuento de Poisson, donde se sabía cuántas veces ocurrió un evento, pero no cuántas veces no ocurrió.
Procesos de modelo que involucran variables de respuesta proporcionales en R especificando un modelo lineal generalizado con familia = binomio. La única complicación es que mientras que con los errores de Poisson podrían simplemente especificar \(family = poisson\), con los errores binomiales se debe dar el número de fracasos así como el número de éxitos en una variable de respuesta de dos vectores. Para hacer esto, une dos vectores usando \(cbind\) en un solo objeto, y, que comprende el número de éxitos y el número de fracasos. El denominador binomial, n, es la muestra total, y:
\[number~of~failures~<-~binomial~denominator~-~number~of~successes\\ y~<-~cbind(number~of~successes, number~of~failures)\]
La forma antigua de modelar este tipo de datos era utilizar el porcentaje de mortalidad como variable de respuesta. Hay cuatro problemas con esto:
Los errores no se distribuyen normalmente.
La varianza no es constante.
La respuesta está limitada (por 1 arriba y 0 abajo).
Al calcular el porcentaje, se pierde información sobre el tamaño de la muestra, n, a partir de la cual se estimó la proporción.
R lleva a cabo una regresión ponderada, utilizando los tamaños de muestra individuales como ponderaciones y la función de enlace logit para garantizar la linealidad. Hay algunos tipos de datos proporcionales, como la cobertura porcentual, que se analizan mejor utilizando modelos lineales convencionales (asumiendo errores normales y varianza constante) después de la transformación de arcoseno. La variable de respuesta, y, medida en radianes,es \[sin^{-1}\sqrt{0.01*p}\], donde \(p\) es el porcentaje de covertura. Sin embargo, si la variable de respuesta toma la forma de un cambio porcentual en alguna medición continua (como el cambio porcentual en el peso al recibir una dieta en particular), entonces, en lugar de transformar los datos en arcoseno, generalmente se trata mejor con:
Análisis de covarianza, utilizando el peso final como variable de respuesta y el peso inicial como covariable.
Especificando la variable de respuesta como una tasa de crecimiento relativa, medida como log (peso final / peso inicial)
Ambos pueden analizarse como modelos lineales con errores normales sin más transformación.
16.1 - Análisis de información en una y dos proporciones
Para comparaciones de una proporción binomial con una constante, use \(binom.test\). Para comparar dos muestras de datos proporcionales, utilice \(prop.test\). Los métodos de este capítulo son necesarios solo para modelos más complejos de datos proporcionales, incluidas las tablas de regresión y contingencia, donde se utilizan \(glm\).
16.2 - Recuento de datos en proporciones
Las transformaciones tradicionales de los datos proporcionales eran arcoseno y probit. La transformación de arcoseno se encargó de la distribución del error, mientras que la transformación probit se utilizó para linealizar la relación entre el porcentaje de mortalidad y la dosis logarítmica en un bioensayo. No hay nada de malo en estas transformaciones y están disponibles dentro de R, pero a menudo es preferible un enfoque más simple, y es probable que produzca un modelo que sea más fácil de interpretar.
La principal dificultad con el modelado de datos proporcionales es que las respuestas están estrictamente limitadas. No hay forma de que el porcentaje de muerte sea mayor al 100% o menor al 0%. Pero si se utilizan técnicas simples como la regresión o el análisis de covarianza, entonces el modelo ajustado podría predecir fácilmente valores negativos o valores superiores al 100%, especialmente si la varianza era alta y muchos de los datos estaban cerca de 0 o cerca del 100%. .
La curva logística se usa comúnmente para describir datos sobre proporciones porque, a diferencia del modelo de línea recta, tiene asíntotas en 0 y 1, por lo que no se pueden predecir proporciones negativas y respuestas de más del 100%. A lo largo de esta discusión se utilizará p para describir la proporción de individuos observados que responden de una manera determinada. Debido a que gran parte de su jerga se deriva de la teoría del juego, los estadísticos llaman a estos éxitos, aunque para un demógrafo que mida las tasas de mortalidad esto puede parecer algo macabro. La proporción de individuos que responden de otras formas (fallas del estadístico) es, por lo tanto, \(1 - p\), y llamará a esta proporción q. La tercera variable es el tamaño de la muestra, n, a partir de la cual se estimó \(p\) (este es el denominador binomial y el número de intentos del estadístico).
Un punto importante sobre la distribución binomial es que la varianza no es constante. De hecho, la varianza de una distribución binomial con media np es:
\[s^2 = npq\]
Entonces la varianza cambia con la varianza de la sioguiente forma:
media <- rnorm(300000)
varianza <- dnorm(media)
ggplot(data.frame(x = media, y = varianza)) +
aes(x = x, y = y) +
geom_point() +
labs(x = "Media", y = "Varianza")
La varianza es baja cuando p es muy alta o muy baja, y la varianza es mayor cuando \(p = q = 0.5\). A medida que p se hace más pequeño, la distribución binomial se acerca cada vez más a la distribución de Poisson. Puede ver por qué esto es así al considerar la fórmula para la varianza del binomio (arriba). Recuerde que para Poisson, la varianza es igual a la media: \(s^{2}= np\). Ahora, a medida que \(p\) se hace más pequeño, \(q\) se acerca cada vez más a 1, por lo que la varianza del binomio converge a la media:
\[s^{2}=npq \approx np ~~~~~~~~~~~~~~~~(q\approx1)\]
16.3 - Posibilidades
El modelo logístico para p en función de x viene dado por:
\[p=\frac{e^{a+bx}}{1+e^{a+bx}}\]
y no hay premios por darse cuenta de que el modelo no es lineal. Pero si x = –∞ entonces p = 0, y si x = + ∞ entonces p = 1, entonces el modelo está estrictamente acotado. Si x = 0, entonces \(p=\frac{e^{a}}{1+e^{a}}\). El truco de linealizar el modelo logístico en realidad implica una transformación muy simple. Es posible que se haya encontrado con la forma en que los corredores de apuestas especifican las probabilidades citando las probabilidades en contra de que un caballo en particular gane una carrera (podrían dar probabilidades de 2 a 1 en un caballo razonablemente bueno o de 25 a 1 en un extraño). Esta es una forma bastante diferente de presentar información sobre probabilidades a la que están acostumbrados los científicos. Por lo tanto, cuando el científico podría establecer una proporción de 0,333 (una probabilidad de ganar en tres), la casa de apuestas daría probabilidades de 2 a 1 (según el conteo de resultados: un éxito contra dos fracasos). En símbolos, esta es la diferencia entre el científico que indica la probabilidad p y la casa de apuestas que indica las probabilidades p / q. Ahora, si se toman las probabilidades p / q y se sustituyen en la fórmula de la logística, se obtiene:
\[\frac{p}{q}=\frac{e^{a+bx}}{1+e^{a+bx}}*[1-\frac{e^{a+bx}}{1+e^{a+bx}}]^{-1}\] el cual se ve terrorifico. Pero con un poco de aritmetica se ve:
\[\frac{p}{q}=\frac{e^{a+bx}}{1+e^{a+bx}}*[1-\frac{e^{a+bx}}{1+e^{a+bx}}]^{-1}=e^{a+bx}\]
Tomar registros naturales y recordar que \(ln(e^x)=x\) simplificará aún más las cosas, de modo que
\[ln(\frac{p}{q})=a+bx\]
Esto da un predictor lineal, \(a + bx\), no para p sino para la transformación logit de p, a saber, \(ln (\frac{p}{q})\). En la jerga de R, el logit es la función de enlace que relaciona el predictor lineal con el valor de p. Aquí se muestran p como función de x (panel izquierdo) y logit (p) como función de x (panel derecho) para el modelo logístico con a = 1 yb = 0,1:
X <- seq(-60,60,by = 0.2)
length(X)## [1] 601P <- pnorm(X, mean = 2.5, sd = 2)
logit <- seq(-6,6,by=0.02)
length(logit)## [1] 601par(mfrow=c(1,2))
plot(x=X,y=logit,ylab="logit = log(p/q)")
plot(x=X,y=P)
En esta etapa, podría preguntarse: “¿por qué no simplemente hacer una regresión lineal de \(ln (\frac{p}{q})\) contra la variable explicativa x?” GLM con errores binomiales tiene tres grandes ventajas aquí:
Permite la varianza binomial no constante.
Se trata del hecho de que los logits para ps cercanos a 0 o 1 son infinitos.
Permite diferencias entre los tamaños de las muestras mediante regresión ponderada.
16.4 Sobredispersión y prueba de hipótesis
Que los errores se distribuyan binomialmente es una suposición, no un hecho. Cuando hay sobredispersión, esta suposición es incorrecta y hay que lidiar con ella.
Todos los diferentes procedimientos estadísticos que se han cumplido en capítulos anteriores también se pueden utilizar con datos sobre proporciones. Se pueden realizar análisis factoriales de varianza, regresión múltiple y una variedad de modelos en los que se ajustan diferentes líneas de regresión en cada uno de varios niveles de uno o más factores. La única diferencia es que evalúan la importancia de los términos sobre la base de chi-cuadrado: el aumento en la desviación escalada que resulta de la eliminación del término del modelo actual.
El punto importante a tener en cuenta es que la prueba de hipótesis con errores binomiales es menos clara que con errores normales. Si bien la aproximación de chi-cuadrado para los cambios en la desviación escalada es razonable para muestras grandes (es decir, más de aproximadamente 30), es más deficiente con muestras pequeñas. Lo más preocupante es el hecho de que se desconoce el grado en que la aproximación es satisfactoria. Esto significa que debe ejercerse un cuidado considerable en la interpretación de las pruebas de hipótesis sobre parámetros, especialmente cuando los parámetros son marginalmente significativos o cuando explican una fracción muy pequeña de la desviación total. Con errores binomiales o de Poisson, no se puede esperar proporcionar valores p exactos para nuestras pruebas de hipótesis.
Cuando se ha obtenido el modelo mínimo adecuado, la desviación escalada residual debe ser aproximadamente igual a los grados de libertad residuales. La sobredispersión ocurre cuando la desviación residual es mayor que los grados de libertad residuales. Hay dos posibilidades: o el modelo está mal especificado o la probabilidad de éxito, \(p\), no es constante dentro de un nivel de tratamiento dado. El efecto de variar p aleatoriamente es aumentar la varianza binomial de \(npq\) a:
\[s^2=npq+n(n-1)\sigma^2\]
Dando lugar a una gran desviación residual. Esto ocurre incluso para modelos que encajarían bien si la variación aleatoria se especificara correctamente.
Una solución simple es asumir que la varianza no es \(npq\) sino \(npq\phi\), donde \(\phi\) es un parámetro de escala desconocido (\(\phi> 1\)). Obtenemos una estimación del parámetro de escala dividiendo el chi-cuadrado de Pearson por los grados de libertad y usamos esta estimación de \(\phi\) para comparar las desviaciones escaladas resultantes. Para lograr esto, se usa family = quasibinomial en lugar de family = binomial cuando hay dispersión excesiva.
Los puntos más importantes a enfatizar en el modelado con errores binomiales son los siguientes:
Cree un objeto de dos columnas para la respuesta, usando cbind para unir los dos vectores que contienen los conteos de éxito y fracaso.
Verifique la sobredispersión (desviación residual mayor que los grados residuales de libertad) y corríjala utilizando family = quasibinomial en lugar de family = binomial si es necesario.
Recuerde que no obtiene valores p exactos con errores binomiales; las aproximaciones de chi-cuadrado son correctas para muestras grandes, pero las muestras pequeñas pueden presentar un problema.
Los valores ajustados son dos conjuntos de recuentos, como la variable de respuesta.
El predictor lineal está en logits (el logaritmo de las probabilidades = \(ln (\frac{p}{q})\)).
Puede volver a transformar de logits (z) a proporciones (p) por p = 1 / [1 + 1 / exp (z)].
16.5 - Aplicaciones
Puede realizar tantos tipos de modelado en un GLM como en un modelo lineal. A continuación se muestran ejemplos de:
Regresión con errores binomiales (variables explicativas continuas).
Análisis de desviación con errores binomiales (variables explicativas categóricas).
Análisis de covarianza con errores binomiales (ambos tipos de variables explicativas).
16.5.1 Regresión logistica con errores binomiales
Este ejemplo se refiere a las proporciones de sexos en los insectos (la proporción de todos los individuos que son machos). En la especie en cuestión, se ha observado que la proporción de sexos es muy variable y se organizó un experimento para ver si la densidad de población estaba involucrada en la determinación de la fracción de machos.
numbers <- read.csv("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/sexratio.txt", header = T,sep = "")
attach(numbers)
head(numbers)## density females males
## 1 1 1 0
## 2 4 3 1
## 3 10 7 3
## 4 22 18 4
## 5 55 22 33
## 6 121 41 80Ciertamente, parece que hay proporcionalmente más machos en alta densidad, pero debería graficar los datos como proporciones para ver esto con mayor claridad:
Densidad <- numbers$density
Hembras <- numbers$females
Machos <- numbers$males
p <- Machos/(Machos+Hembras)
log_den <- log(Densidad)
dat_ins <- data.frame(Densidad,Hembras,Machos,p,log_den)
fig_den <- plot_ly(data = dat_ins,x = Densidad,y = p)
fig_den <- fig_den %>% layout(
xaxis = list(title = "Densidad"),
yaxis = list(title = "Proporción de Machos"))
fig_den## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-modefig_den_1 <- plot_ly(data = dat_ins,x = log_den,y = p)
fig_den_1 <- fig_den_1 %>% layout(
xaxis = list(title = "log(Densidad)"),
yaxis = list(title = "Proporción de Machos"))
fig_den_1## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plot.ly/r/reference/#scatter## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plot.ly/r/reference/#scatter-modeEvidentemente, es probable que una transformación logarítmica de la variable explicativa mejore el ajuste del modelo. lo verá en un momento.
La pregunta es si el aumento de la densidad de población conduce a un aumento significativo en la proporción de machos en la población o, más brevemente, si la proporción de sexos depende de la densidad. Ciertamente, se ve en la trama como si lo fuera.
La variable de respuesta es un par de recuentos coincidentes que deseamos analizar como datos proporcionales usando un GLM con errores binomiales. Primero, usamos cbind para unir los vectores de conteos masculinos y femeninos en un solo objeto que será la respuesta en nuestro análisis:
y <- cbind(Machos,Hembras)Esto significa que \(y\) se interpretará en el modelo como la proporción de todos los individuos que eran hombres. El modelo se especifica así:
model_16 <- glm(y~Densidad,binomial)Esto dice que el objeto llamado modelo obtiene un modelo lineal generalizado en el que \(y\) (la proporción de sexos) se modela como una función de una única variable explicativa continua (llamada densidad), utilizando una distribución de error de la familia binomial. La salida se ve así:
summary(model_16)##
## Call:
## glm(formula = y ~ Densidad, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.4619 -1.2760 -0.9911 0.5742 1.8795
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0807368 0.1550376 0.521 0.603
## Densidad 0.0035101 0.0005116 6.862 6.81e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 71.159 on 7 degrees of freedom
## Residual deviance: 22.091 on 6 degrees of freedom
## AIC: 54.618
##
## Number of Fisher Scoring iterations: 4La tabla del modelo tiene el mismo aspecto que para una regresión sencilla. El primer parámetro en la tabla de Coeficientes es la intersección y el segundo es la pendiente del gráfico de la proporción de sexos contra la densidad de población. La pendiente es significativamente más pronunciada que cero (proporcionalmente más machos con mayor densidad de población; \(p = 6,81 × 10-12\)), pero hay una sobredispersión sustancial (la desviación residual = 22,091 es mucho mayor que la d.f. residual = 6). Podemos ver si la transformación logarítmica de la variable explicativa mejora esto:
model17 <- glm(y~log_den,binomial)
summary(model17)##
## Call:
## glm(formula = y ~ log_den, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9697 -0.3411 0.1499 0.4019 1.0372
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.65927 0.48758 -5.454 4.92e-08 ***
## log_den 0.69410 0.09056 7.665 1.80e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 71.1593 on 7 degrees of freedom
## Residual deviance: 5.6739 on 6 degrees of freedom
## AIC: 38.201
##
## Number of Fisher Scoring iterations: 4Esta es una gran mejora, por lo que se adoptará. En el modelo con log (densidad) no hay evidencia de sobredispersión (desviación residual = 5.67 en 6 d.f.), mientras que la falta de ajuste introducida por la curvatura en el primer modelo causó una sobredispersión sustancial (desviación residual = 22.09 en 6 d.f.).
La verificación del modelo implica el uso de un graficador para el modelo 17. Como verá, no hay un patrón en los residuales contra los valores ajustados y la gráfica normal es razonablemente lineal. Punto no. 4 es muy influyente (tiene una gran distancia de Cook), pero el modelo sigue siendo significativo con este punto omitido.
Se concluye que la proporción de insectos que son machos aumenta significativamente al aumentar la densidad, y que el modelo logístico se linealiza por transformación logarítmica de la variable explicativa (densidad de población). Termina dibujando la línea ajustada a través del diagrama de dispersión:
xv_4 <- seq(0,6,0.01)
yv_4 <- predict(model17,list(exp(xv_4)),type = "response")
plot(log_den,p)
#lines(xv_4,yv_4)
#df_gg <- data.frame(xv_4,yv_5)
#ggplot(dat_ins, aes(x = log_den, y = p)) +
# geom_point() +
#geom_line(color ='red',data=df_gg, aes(x=xv_4,y=yv_5))+
#labs(x="log(Densidad",y="Proporción machos")Nótese el uso de type = “response” para retrotransformar de la escala logit a la escala de proporción en forma de S.
16.5.2 - Estimando LD50 y LD90 desde la información de biomasa
Los datos consisten en el número de muertos y el tamaño del lote inicial para cinco dosis de aplicación de plaguicidas, y se desea saber qué dosis mata al 50% de los individuos (o al 90% o 95%, según sea necesario). El problema estadístico complicado es que se está usando un valor de \(y\) (50% muerto) para predecir un valor de \(x\) (la dosis relevante) y calcular un error estándar en el eje x.
data_1 <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/bioassay.txt",header = T)
attach(data_1)
names(data_1)## [1] "dose" "dead" "batch"La regresión logistica es llevada a cabo en el camino usual
y_1 <- cbind(dead,batch-dead)
model18 <- glm(y_1~log(dose),binomial)Entoces la función \(dose.p\) de la librería MASS se corre con el objeto del modelo, especificando las proporciones “matadas” para los cuales se quiere el \(log(doses)\) predicho (\(p\)=0.5 es el predeterminado para LD50)
dose.p(model18,p=c(0.5,0.9,0.95))## Dose SE
## p = 0.50: 2.306981 0.07772065
## p = 0.90: 3.425506 0.12362080
## p = 0.95: 3.805885 0.15150043A pesar de la etiqueta “Dosis”, la salida muestra los registros de las dosis asociadas con muertes de LD50, LD90 y LD95, junto con sus errores estándar.
16.5.3 - Proporción de información con variables explicativas categóricas
El siguiente ejemplo se refiere a la germinación de semillas de dos genotipos de la planta parásita Orobanche y dos extractos de plantas hospedantes (frijol y pepino) que se utilizaron para estimular la germinación. Es un análisis factorial bidireccional de la desviación.
germination <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/germination.txt")
attach(germination)
names(germination)## [1] "count" "sample" "Orobanche" "extract"count_1 <- germination$countEl recuento es el número de semillas que germinaron de un lote de \(tamaño = muestra\). Entonces, el número que no germinó es \(muestra - recuento\), y se construye el vector de respuesta así:
y_2 <- cbind(count_1,sample-count_1)Cada uno de las variables explicativas categóricas tiene 2 niveles
Orobanche <- as.factor(Orobanche)
levels(Orobanche)## [1] "a73" "a75"extract <- as.factor(extract)
levels(extract)## [1] "bean" "cucumber"Se quiere contrastar la hipótesis de que no existe interacción entre el genotipo Orobanche (a73 o a75) y el extracto vegetal (frijol o pepino) en la tasa de germinación de las semillas. Esto requiere un análisis factorial usando el operador asterisco \(*\) como este:
model_19 <- glm(y_2~Orobanche*extract,binomial)
summary(model_19)##
## Call:
## glm(formula = y_2 ~ Orobanche * extract, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.01617 -1.24398 0.05995 0.84695 2.12123
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4122 0.1842 -2.238 0.0252 *
## Orobanchea75 -0.1459 0.2232 -0.654 0.5132
## extractcucumber 0.5401 0.2498 2.162 0.0306 *
## Orobanchea75:extractcucumber 0.7781 0.3064 2.539 0.0111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 98.719 on 20 degrees of freedom
## Residual deviance: 33.278 on 17 degrees of freedom
## AIC: 117.87
##
## Number of Fisher Scoring iterations: 4A primera vista, parece que hay una interacción muy significativa (\(p = 0,0111\)). Pero es necesario comprobar que el modelo es sólido. Lo primero que hay que comprobar es la sobredispersión. La desviación residual es 33.278 en 17 d.f., por lo que el modelo está bastante sobredispersado:
33.278/17## [1] 1.957529El factor de sobredispersión es casi 2. La forma más sencilla de tener esto en cuenta es utilizar lo que se llama un “parámetro de escala empírica” para reflejar el hecho de que los errores no son binomiales como se suponía, sino que eran mayores que esto (es decir, dispersos en exceso) en un factor de 1,9576. Se reajusta el modelo usando errores cuasi-binomiales para dar cuenta de la sobredispersión:
model_20 <- glm(y_2~Orobanche*extract,quasibinomial)
model_21 <- glm(y_2~Orobanche+extract,quasibinomial)La única diferencia es que se usa una prueba F en lugar de una prueba de chi-cuadrado para comparar los modelos original y simplificado porque ahora se han estimado dos parámetros del modelo (la media más el parámetro de escala empírica):
anova(model_20,model_21,test = "F")## Analysis of Deviance Table
##
## Model 1: y_2 ~ Orobanche * extract
## Model 2: y_2 ~ Orobanche + extract
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 17 33.278
## 2 18 39.686 -1 -6.4081 3.4418 0.08099 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ahora se observa que la interacción no es significativa (\(p=0.081\)). No hay evidencia convicente de que los diferentes genotipos de Orobanche, responda diferente ante los dos extractos de las plantas.
El siguiente paso es ver si es posible una mayor simplificación del modelo:
anova(model_21,test="F")## Analysis of Deviance Table
##
## Model: quasibinomial, link: logit
##
## Response: y_2
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 20 98.719
## Orobanche 1 2.544 19 96.175 1.1954 0.2887
## extract 1 56.489 18 39.686 26.5412 6.692e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model_22 <- glm(y_2~extract,quasibinomial)
anova(model_21,model_22,test="F")## Analysis of Deviance Table
##
## Model 1: y_2 ~ Orobanche + extract
## Model 2: y_2 ~ extract
## Resid. Df Resid. Dev Df Deviance F Pr(>F)
## 1 18 39.686
## 2 19 42.751 -1 -3.065 1.4401 0.2457No hay justificación suficiente para retener a Orobanche dentro del modelo. Entonces la minima adecuación del modelo contiene solo dos parámetros:
coef(model_22)## (Intercept) extractcucumber
## -0.5121761 1.0574031¿Qué significan exactamente estos dos números? Recuerde que los coeficientes son del predictor lineal. Están en la escala transformada, así que debido a que estamos usando errores cuasi-binomiales, están en logits (\(ln (\frac{p}{1 - p})\)). Para convertirlos en las tasas de germinación para los dos extractos de plantas se requiere un pequeño cálculo. Para pasar de un logit \(x\) a una proporción \(p\), es necesario calcular
\[p=\frac{1}{1+1/e^x}\]
Entonces el primer valor de \(x\) es -0.5122, por lo que se calcula
1/(1+1/exp(-0.5122))## [1] 0.3746779Esto dice que la tasa de germinación media de las semillas con el primer extracto de planta (frijol) fue del 37%. ¿Qué pasa con el parámetro de extracción (1.0574)? Recuerde que con las variables explicativas categóricas los valores de los parámetros son diferencias entre medias. Entonces, para obtener la segunda tasa de germinación, agregue 1.057 a la intersección antes de la transformación inversa:
1/(1+1/exp(-0.5122+1.0574))## [1] 0.6330212Esto dice que la tasa de germinación fue casi el doble (63%) con el segundo extracto de planta (pepino). Evidentemente se quiere generalizar este proceso, y también agilizar los cálculos de las proporciones medias estimadas. Puede usar \(predict\) para ayudar aquí, porque \(type = "response"\) hace predicciones en la escala transformada hacia atrás automáticamente:
tapply(predict(model_22,type = "response"),extract,mean)## bean cucumber
## 0.3746835 0.6330275Es interesante comparar estas cifras con los promedios de las proporciones brutas. Primero es necesario calcular la proporción de germinación, p, en cada muestra:
p_1 <- count_1/sampleLuego se pueden encontrar las tasas de germinación promedio de cada extracto:
tapply(p_1,extract,mean)## bean cucumber
## 0.3487189 0.6031824Ve que esto da diferentes respuestas. No muy diferente en este caso, pero diferente de todos modos. La forma correcta de promediar los datos de proporción es sumar los recuentos totales para los diferentes niveles de resumen, y solo entonces convertirlos en proporciones:
tapply(count_1,extract,sum)## bean cucumber
## 148 276Esto significa que 148 semillas germinaron con extracto de frijol y 276 con pepino. Pero, ¿cuántas semillas estuvieron involucradas en cada caso?
tapply(sample,extract,sum)## bean cucumber
## 395 436Esto significa que 395 semillas fueron tratadas con extracto de frijol y 436 semillas fueron tratadas con pepino. Entonces, las respuestas que queremos son 148/395 y 276/436 (es decir, las proporciones medias correctas). Podemos automatizar el cálculo de esta manera:
as.vector(tapply(count_1,extract,sum))/as.vector(tapply(sample,extract,sum))## [1] 0.3746835 0.6330275Estas son las proporciones medias correctas que fueron producidas por el GLM. La moraleja aquí es que se calcula el promedio de proporciones usando recuentos totales y muestras totales y no promediando las proporciones brutas.
Promediando Proporciones
aquí un ejemplo de lo que NO hay que hacer. Se tienen 4 proporciones
\[0.2,~0.17,~0.2,~0.53\]
Así que seguramente es solo sumarlos y dividirlos por 4. Esto da \(1.1 / 4 = 0.275\). ¡Incorrecto! Y no por un poquito. Es necesario observar los recuentos en los que se basaron las proporciones. Estos resultan ser:
\[1/5,~1/6,~2/10,~53/100\]
La forma correcta de promediar las proporciones es sumar el recuento total de éxitos (\(1 + 1 + 2 + 53 = 57\)) y dividirlo por el número total de muestras (\(5 + 6 + 10 + 100 = 121\)). La proporción media correcta es \(57/121 = 0,4711\). Esto es casi el doble de nuestra respuesta incorrecta (arriba).
16.7 - Resumen del modelado con recuento de datos de proporciones
Cree un vector de respuesta de dos columnas que contenga los éxitos y los fracasos.
Use \(glm\) con \(family = binomial\) (puede omitir \(family =\)).
Ajuste el modelo máximo.
Prueba de sobredispersión.
Si encuentra una dispersión excesiva, utilice errores cuasibinomiales en lugar de binomiales.
Comience la simplificación del modelo eliminando los términos de interacción.
Eliminar efectos principales no significativos
Utilice plot para obtener sus diagnósticos de verificación de modelo.
Realice una transformación inversa usando \(predict\) con la opción \(type = "response"\) para obtener medias.
Análisis de covarianza con información binomial
Pasemos ahora a un ejemplo sobre la floración en cinco variedades de plantas perennes. Los individuos replicados en un diseño completamente aleatorizado se rociaron con una de seis dosis de una mezcla controlada de promotores del crecimiento. Después de 6 semanas, las plantas se calificaron como con o sin floración. El recuento de individuos en floración constituye la variable de respuesta. Este es un ANCOVA porque tiene variables explicativas tanto continuas (dosis) como categóricas (variedad). Se usa regresión logística porque la variable de respuesta es un recuento (flores) que se puede expresar como una proporción (flores / número).
props <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/flowering.txt",header = T)
attach(props)
names(props)## [1] "flowered" "number" "dose" "variety"y_3 <- cbind(flowered,number-flowered)
pf <- flowered/number
pfc <- split(pf,variety)
dc <- split(dose,variety)
dc## $A
## [1] 1 4 8 16 32 0
##
## $B
## [1] 1 4 8 16 32 0
##
## $C
## [1] 1 4 8 16 32 0
##
## $D
## [1] 1 4 8 16 32 0
##
## $E
## [1] 1 4 8 16 32 0pfc## $A
## [1] 0.0000000 0.0000000 0.4000000 0.8181818 1.0000000 0.0000000
##
## $B
## [1] 0.00000000 0.20000000 0.50000000 0.90000000 0.50000000 0.08333333
##
## $C
## [1] 0.14285714 0.06666667 0.17647059 0.25000000 1.00000000 0.00000000
##
## $D
## [1] 0.11111111 0.15789474 0.53571429 0.73076923 0.77777778 0.06666667
##
## $E
## [1] 0.0000000 0.0000000 0.1578947 0.7500000 1.0000000 0.0000000plot(dose,pf,type="n",ylab="Proportion flowered")
points(jitter(dc[[1]]),jitter(pfc[[1]]),pch=21,col="blue",bg="red")
points(jitter(dc[[2]]),jitter(pfc[[2]]),pch=21,col="blue",bg="green")
points(jitter(dc[[3]]),jitter(pfc[[3]]),pch=21,col="blue",bg="yellow")
points(jitter(dc[[4]]),jitter(pfc[[4]]),pch=21,col="blue",bg="green3")
points(jitter(dc[[5]]),jitter(pfc[[5]]),pch=21,col="blue",bg="brown")
## Cambiar gráficoTenga en cuenta el uso de split para separar las diferentes variedades, de modo que podamos trazarlas con diferentes símbolos, y de jitter para evitar que los valores repetidos se oculten entre sí. Es evidente que existe una diferencia sustancial entre las variedades de plantas en su respuesta al estimulante de la floración. El modelado procede de la forma habitual. Comenzamos ajustando el modelo máximo con diferentes pendientes e intersecciones para cada variedad (estimando diez parámetros en total):
model_23 <- glm(y_3~dose*variety,binomial)
summary(model_23)##
## Call:
## glm(formula = y_3 ~ dose * variety, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.6648 -1.1200 -0.3769 0.5735 3.3299
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.59165 1.03215 -4.449 8.64e-06 ***
## dose 0.41262 0.10033 4.113 3.91e-05 ***
## varietyB 3.06197 1.09317 2.801 0.005094 **
## varietyC 1.23248 1.18812 1.037 0.299576
## varietyD 3.17506 1.07516 2.953 0.003146 **
## varietyE -0.71466 1.54849 -0.462 0.644426
## dose:varietyB -0.34282 0.10239 -3.348 0.000813 ***
## dose:varietyC -0.23039 0.10698 -2.154 0.031274 *
## dose:varietyD -0.30481 0.10257 -2.972 0.002961 **
## dose:varietyE -0.00649 0.13292 -0.049 0.961057
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 303.350 on 29 degrees of freedom
## Residual deviance: 51.083 on 20 degrees of freedom
## AIC: 123.55
##
## Number of Fisher Scoring iterations: 5El modelo exhibe una sobredispersión sustancial, pero esto probablemente se deba a una mala selección del modelo en lugar de una variabilidad adicional no medida. Investiguemos esto trazando las curvas ajustadas a través del diagrama de dispersión.
# gráfico p.p. 664Como puede ver, el modelo es razonable para dos de los genotipos (A y E, representados por líneas rojas y marrones, respectivamente), moderado para un genotipo (C, amarillo) pero muy pobre para dos de ellos, B (inferior, línea verde claro) y D (línea superior, verde oscuro). Para ambos últimos, el modelo sobreestima en gran medida la proporción de floración a dosis cero, y para el genotipo B parece haber cierta inhibición de la floración a la dosis más alta porque el gráfico cae del 90% de floración a la dosis 16 a solo el 50% a la dosis. 32. La variedad D parece ser asintótica con menos del 100% de floración.
tapply(pf,list(dose,variety),mean)## A B C D E
## 0 0.0000000 0.08333333 0.00000000 0.06666667 0.0000000
## 1 0.0000000 0.00000000 0.14285714 0.11111111 0.0000000
## 4 0.0000000 0.20000000 0.06666667 0.15789474 0.0000000
## 8 0.4000000 0.50000000 0.17647059 0.53571429 0.1578947
## 16 0.8181818 0.90000000 0.25000000 0.73076923 0.7500000
## 32 1.0000000 0.50000000 1.00000000 0.77777778 1.0000000Estos fallos del modelo deberían centrar la atención para el trabajo futuro. La moraleja es que el hecho de que tengamos datos proporcionales no significa que los datos necesariamente estarán bien descritos por el modelo logístico. Por ejemplo, para describir la respuesta del genotipo B, el modelo debería tener una joroba, en lugar de una asíntota en p = 1 para dosis grandes. Es esencial observar de cerca los datos, tanto con gráficos como con tablas, antes de aceptar el resultado del modelo. La elección del modelo es muy importante. La logística fue una mala elección para dos de las cinco variedades en este caso.
16.9 - Convirtiendo tablas de contingencia complejas a proporciones
En esta sección se muestra cómo eliminar la necesidad de todas las variables molestas que están involucradas en el modelado de tablas de contingencia complejas (los totales de filas y columnas que deben incluirse en todos los modelos). Puede hacer esto cuando la respuesta puede reformularse como una proporción binaria (una opción de entre dos contingencias). Por ejemplo, en el caso de las lagartijas de Schoener, que resultó tan difícil de analizar como una tabla de contingencia compleja (ver p. 610), se puede trabajar con la proporción de todas las lagartijas que son Anolis grahamii como variable de respuesta, en lugar de analizar los recuentos de los números de A. grahamii y A. opalinus por separado. Esto tiene la gran ventaja de no requerir que ninguna de las variables molestas se incluya en el modelo. Sin embargo, la técnica no funcionaría si tuviéramos tres especies de lagartos. Entonces tendría que ceñirse al complejo modelo de tablas de contingencia.
lizards <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/lizards.txt",header = T)
attach(lizards)
names(lizards)## [1] "n" "sun" "height" "perch" "time" "species"head(lizards)## n sun height perch time species
## 1 20 Shade High Broad Morning opalinus
## 2 13 Shade Low Broad Morning opalinus
## 3 8 Shade High Narrow Morning opalinus
## 4 6 Shade Low Narrow Morning opalinus
## 5 34 Sun High Broad Morning opalinus
## 6 31 Sun Low Broad Morning opalinusPrimero, es necesario estar absolutamente seguro de que todas las variables explicativas estén exactamente en el mismo orden para ambas especies de lagartos. La razón de esto es que se van a unir los recuentos de una de las especies de lagartos en la mitad del marco de datos que contiene los recuentos de otras especies y todas las variables explicativas. Cualquier error aquí sería desastroso porque el recuento estaría alineado con la combinación incorrecta de variables explicativas, y el análisis sería incorrecto y carecería de sentido.
species <- as.factor(lizards$species)
sun <- as.factor(lizards$sun)
height <- as.factor(lizards$height)
perch <- as.factor(lizards$perch)
time <- as.factor(lizards$time)
sorted <- lizards[order(species,sun,height,perch,time),]
head(sorted)## n sun height perch time species
## 41 4 Shade High Broad Afternoon grahamii
## 33 1 Shade High Broad Mid.day grahamii
## 25 2 Shade High Broad Morning grahamii
## 43 3 Shade High Narrow Afternoon grahamii
## 35 1 Shade High Narrow Mid.day grahamii
## 27 3 Shade High Narrow Morning grahamiiA continuación, es necesario extraer la mitad superior de este marco de datos.
short <- sorted[1:24,]Tenga en cuenta que esta operación ha perdido todos los datos de A. opalinus. Además, el nombre de la variable de la izquierda, n, ya no es apropiado. Es el recuento de \(A. grahamii\), por lo que deberíamos renombrarlo Ag, digamos (con la intención de agregar otra columna llamada Ao a su debido tiempo para contener los recuentos de \(A. opalinus\)):
names(short)[1] <- "Ag"
names(short)## [1] "Ag" "sun" "height" "perch" "time" "species"La variable de la derecha, especie, es redundante ahora (todas las entradas son \(A. grahamii\)), por lo que debería eliminarla:
short <- short[,-6]
head(short)## Ag sun height perch time
## 41 4 Shade High Broad Afternoon
## 33 1 Shade High Broad Mid.day
## 25 2 Shade High Broad Morning
## 43 3 Shade High Narrow Afternoon
## 35 1 Shade High Narrow Mid.day
## 27 3 Shade High Narrow MorningLos recuentos de cada fila de A. opalinus están en la variable llamada n en la mitad inferior del marco de datos llamado \(sorted\). Los extraemos así:
sorted$n[25:48]## [1] 4 8 20 5 4 8 12 8 13 1 0 6 18 69 34 8 60 17 13 55 31 4 21 12La idea es crear un nuevo marco de datos con estos recuentos de \(A. opalinus\) alineados junto con los recuentos coincidentes de \(A.grahamii\):
new.lizards <- data.frame(sorted$n[25:48],short)La primera variable necesita un nombre informativo, como \(Ao\):
names(new.lizards)[1] <- "Ao"
head(new.lizards)## Ao Ag sun height perch time
## 41 4 4 Shade High Broad Afternoon
## 33 8 1 Shade High Broad Mid.day
## 25 20 2 Shade High Broad Morning
## 43 5 3 Shade High Narrow Afternoon
## 35 4 1 Shade High Narrow Mid.day
## 27 8 3 Shade High Narrow MorningEso completa la edición del marco de datos. Tenga en cuenta, sin embargo, que tenemos tres marcos de datos, todas de configuraciones diferentes, pero cada una contiene los mismos nombres de variable (sol, altura, posición y tiempo): observe los objetos () y busque (). Necesitamos hacer algunas tareas de limpieza:
lizards <- as.factor(lizards)
df_lizards <- as.data.frame(lizards)
#detach(df_lizards)
rm(short,sorted)
attach(new.lizards)16.9.1 - Analizando lagartijas de Schoener como una información de proporción
Con los preliminares anteriores, aquí están los nombres de las variables:
names(new.lizards)## [1] "Ao" "Ag" "sun" "height" "perch" "time"La variable respuesta es un obejto de dos columnas el cual contiene el recuento de las dos especies.
y_4 <- cbind(Ao,Ag)Empiece por ajustar el modelo saturado que contiene todas las interacciones posibles:
sun <- as.vector(sun)
height <- as.vector(height)
time <- as.vector(time)
perch <- as.factor(perch)
y_4 <- as.factor(y_4)
model_24 <- glm(y_4 ~ (time*sun*height*perch),binomial)Como no hay variables molestas, podemos usar \(step\) directamente para comenzar la simplificación del modelo:
model_25 <- step(model_24)## Start: AIC=59.09
## y_4 ~ (time * sun * height * perch)
##
## Df Deviance AIC
## - time:sun:height:perch 2 11.09 55.09
## <none> 11.09 59.09
##
## Step: AIC=55.09
## y_4 ~ time + sun + height + perch + time:sun + time:height +
## sun:height + time:perch + sun:perch + height:perch + time:sun:height +
## time:sun:perch + time:height:perch + sun:height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - time:height:perch 2 11.09 51.09
## - time:sun:height 2 11.09 51.09
## - time:sun:perch 2 11.09 51.09
## - sun:height:perch 1 11.09 53.09
## <none> 11.09 55.09
##
## Step: AIC=51.09
## y_4 ~ time + sun + height + perch + time:sun + time:height +
## sun:height + time:perch + sun:perch + height:perch + time:sun:height +
## time:sun:perch + sun:height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - time:sun:height 2 11.09 47.09
## - time:sun:perch 2 11.09 47.09
## - sun:height:perch 1 11.09 49.09
## <none> 11.09 51.09## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=47.09
## y_4 ~ time + sun + height + perch + time:sun + time:height +
## sun:height + time:perch + sun:perch + height:perch + time:sun:perch +
## sun:height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - time:height 2 11.09 43.09
## - time:sun:perch 2 11.09 43.09
## - sun:height:perch 1 11.09 45.09
## <none> 11.09 47.09## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=43.09
## y_4 ~ time + sun + height + perch + time:sun + sun:height + time:perch +
## sun:perch + height:perch + time:sun:perch + sun:height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - time:sun:perch 2 11.09 39.09
## - sun:height:perch 1 11.09 41.09
## <none> 11.09 43.09## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=39.09
## y_4 ~ time + sun + height + perch + time:sun + sun:height + time:perch +
## sun:perch + height:perch + sun:height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - time:perch 2 11.090 35.090
## - sun:height:perch 1 11.090 37.090
## <none> 11.090 39.090
## - time:sun 2 21.504 45.504## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=35.09
## y_4 ~ time + sun + height + perch + time:sun + sun:height + sun:perch +
## height:perch + sun:height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - sun:height:perch 1 11.090 33.090
## <none> 11.090 35.090
## - time:sun 2 22.523 42.523## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=33.09
## y_4 ~ time + sun + height + perch + time:sun + sun:height + sun:perch +
## height:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - height:perch 1 11.090 31.090
## - sun:perch 1 11.090 31.090
## - sun:height 1 12.234 32.234
## <none> 11.090 33.090
## - time:sun 2 22.523 40.523## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=31.09
## y_4 ~ time + sun + height + perch + time:sun + sun:height + sun:perch## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred## Df Deviance AIC
## - sun:perch 1 11.090 29.090
## - sun:height 1 12.873 30.873
## <none> 11.090 31.090
## - time:sun 2 23.068 39.068## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred##
## Step: AIC=29.09
## y_4 ~ time + sun + height + perch + time:sun + sun:height
##
## Df Deviance AIC
## <none> 11.090 29.090
## - sun:height 1 13.388 29.388
## - perch 1 17.995 33.995
## - time:sun 2 23.769 37.769Se ha visto que AIC es muy generoso al dejar términos en el modelo que eliminaríamos sin piedad. Para empezar, necesitamos probar si hubiéramos mantenido las dos interacciones de tres vías y las cinco interacciones bidireccionales:
model_26 <- update(model_25,~.-height:perch:time)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredmodel_27 <- update(model_25,~.-sun:height:perch)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredanova(model_25,model_26,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ time + sun + height + perch + time:sun + sun:height
## Model 2: y_4 ~ time + sun + height + perch + time:sun + sun:height
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 39 11.09
## 2 39 11.09 0 0fue cerca pero no es significante
anova(model_25,model_27,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ time + sun + height + perch + time:sun + sun:height
## Model 2: y_4 ~ time + sun + height + perch + time:sun + sun:height
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 39 11.09
## 2 39 11.09 0 0No, no mantendría ninguna de esas interacciones de tres vías. ¿Qué pasa con las interacciones bidireccionales? es necesario comenzar con un modelo base más simple que el modelo 25:
model_28 <- glm(y_4~(sun+height+perch+time)^2-sun:time,binomial)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredDebe eliminar cada una de las interacciones bidireccionales por separado, comparándolas con el modelo_28 que contiene todas las interacciones bidireccionales:
model_29 <- update(model_28,~.-sun:height)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredanova(model_28,model_29,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y_4 ~ sun + height + perch + time + sun:perch + height:perch +
## height:time + perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 35 21.232
## 2 36 21.504 -1 -0.27267 0.6015model_30 <- update(model_28,~.-sun:perch)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredanova(model_28,model_30,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y_4 ~ sun + height + perch + time + sun:height + height:perch +
## height:time + perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 35 21.232
## 2 36 21.870 -1 -0.63856 0.4242model_31 <- update(model_28,~.-height:perch)## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredanova(model_28,model_31,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y_4 ~ sun + height + perch + time + sun:height + sun:perch +
## height:time + perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 35 21.232
## 2 36 21.658 -1 -0.42617 0.5139model_32 <- update(model_28,~.-time:perch)
anova(model_28,model_32,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y_4 ~ sun + height + perch + time + sun:height + sun:perch +
## height:perch + height:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 35 21.232
## 2 37 21.870 -2 -0.63856 0.7267model_33 <- update(model_28,~.-time:height)
anova(model_28,model_33,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ (sun + height + perch + time)^2 - sun:time
## Model 2: y_4 ~ sun + height + perch + time + sun:height + sun:perch +
## height:perch + perch:time
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 35 21.232
## 2 37 21.504 -2 -0.27267 0.8726Por lo tanto, no necesitamos ninguna de las interacciones bidireccionales. ¿Qué pasa con los efectos principales?
model_34 <- glm(y_4~sun+height+perch+time,binomial)
summary(model_34)##
## Call:
## glm(formula = y_4 ~ sun + height + perch + time, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.25922 0.00005 0.19006 0.39536 1.00916
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 20.138 2449.160 0.008 0.9934
## sunSun 1.031 1.056 0.976 0.3291
## heightLow 1.031 1.056 0.976 0.3291
## perchNarrow -2.113 1.230 -1.718 0.0858 .
## timeMid.day -18.646 2449.160 -0.008 0.9939
## timeMorning -17.615 2449.160 -0.007 0.9943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 36.170 on 47 degrees of freedom
## Residual deviance: 24.449 on 42 degrees of freedom
## AIC: 36.449
##
## Number of Fisher Scoring iterations: 18Todos los efectos principales son importantes y, por tanto, deben mantenerse. Solo un último punto: es posible que no necesitemos los tres niveles para el tiempo, ya que el resumen sugiere que el mediodía y la mañana no son significativamente diferentes (diferencia del parámetro de 0,9639 - 0,7368 = 0,2271, con un error estándar de la diferencia de 0,29). Los agrupamos en un nuevo factor llamado t2:
t3 <- time
time <- as.factor(time)
levels(t3)[c(2)] <- "other"
levels(t3)[c(3)] <- "other1"
levels(t3)## [1] NA "other" "other1"model_35 <- glm(y_4~sun+height+perch+t3,binomial)
anova(model_34,model_35,test = "Chi")## Analysis of Deviance Table
##
## Model 1: y_4 ~ sun + height + perch + time
## Model 2: y_4 ~ sun + height + perch + t3
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 42 24.449
## 2 42 24.449 0 0Un modelo con solo dos horas del día no es significativamente peor que un modelo con tres.
summary(model_35)##
## Call:
## glm(formula = y_4 ~ sun + height + perch + t3, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.25922 0.00005 0.19006 0.39536 1.00916
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 20.138 2449.160 0.008 0.9934
## sunSun 1.031 1.056 0.976 0.3291
## heightLow 1.031 1.056 0.976 0.3291
## perchNarrow -2.113 1.230 -1.718 0.0858 .
## t3Mid.day -18.646 2449.160 -0.008 0.9939
## t3Morning -17.615 2449.160 -0.007 0.9943
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 36.170 on 47 degrees of freedom
## Residual deviance: 24.449 on 42 degrees of freedom
## AIC: 36.449
##
## Number of Fisher Scoring iterations: 18Todos los parámetros son significativos, por lo que este es el modelo mínimo adecuado. No hay evidencia de sobredispersión. Solo hay cinco parámetros y el modelo no contiene variables molestas (compárelo con el modelo de tabla de contingencia masiva summary.lm (modelo1)). La interpretación ecológica es sencilla: las dos especies de lagartos difieren significativamente en sus nichos en todos los ejes de nichos que se midieron. Sin embargo, no hubo interacciones significativas (no sucedía nada sutil, como cambiar el tamaño de las perchas en diferentes momentos del día).
Capítulo 26 - Spatial statistics - Estadística espacial (pp 835)
Hay tres tipos de problemas que puede abordar con las estadísticas espaciales:
procesos puntuales (ubicaciones y patrones espaciales de individuos)
Mapas de una variable respuesta continua (kriging)
Respuestas espaciales explicitas afectadas por la indetidad, tamaño y proximidad de los vecinos.
26.1 Procesos puntuales
Hay tres clases amplias de patrón espacial en un continuo desde la regularidad completa (hexágonos espaciados uniformemente donde cada individuo está a la misma distancia de su vecino más cercano) hasta la agregación completa (todos los individuos agrupados en un solo grupo): los llamamos regulares, aleatorios y patrones agregados y se ven así:
regular_x <- c(2,2,3,3,4,4,5,5,6,6,7,7)
regular_y <- c(3,7,2,5,3,7,2,5,3,7,2,5)
random_x <- runif(14,min = 1,max = 7)
random_y <- runif(14,min = 1,max = 7)
agg_x <- c(2,2.5,2.5,2.5,3,4,4.5,6,6.5,6.5,6.5,6.7,7,7,7)
agg_y <- c(2,1,2,3,2.5,4.5,4,4.5,4,4.7,5,5,4,4.7,5)
par(mfrow=c(1,3))
plot(regular_x,regular_y,ylim = c(0,8),xlim = c(1,8),xlab = "",ylab = "",cex=2,pch=16, main = "Regular")
plot(random_x,random_y,ylim = c(0,8),xlim = c(1,8),xlab = "",ylab = "",cex=2,pch=16, main = "Random")
plot(agg_x,agg_y,ylim = c(0,8),xlim = c(1,8),xlab = "",ylab = "",cex=2,pch=16, main = "aggregated")
En su forma más simple, los datos consisten en conjuntos de coordenadas \(x\) y \(y\) dentro de algún marco de muestreo, como un cuadrado o un círculo, en el que se han mapeado los individuos. La primera pregunta es a menudo si existe alguna evidencia que permita el rechazo de la hipótesis nula de aleatoriedad espacial completa (RSC). En un patrón aleatorio, la distribución de cada individuo es completamente independiente de la distribución de todos los demás. Los individuos no se inhiben ni se promueven entre sí. En un patrón regular, los individuos están más espaciados que en uno aleatorio, presumiblemente debido a algún mecanismo (como la competencia) que elimina a los individuos que están demasiado juntos. En un patrón agregado, los individuos están más agrupados que en uno aleatorio, presumiblemente debido a algún proceso como la reproducción con dispersión limitada, o debido a la heterogeneidad espacial subyacente (por ejemplo, parches buenos y parches malos).
Los recuentos de individuos dentro de áreas de muestra (cuadrantes) se pueden analizar comparando la distribución de frecuencia de los recuentos con una distribución de Poisson con la misma media. Los patrones espaciales agregados (en los que la varianza es mayor que la media) a menudo se describen bien mediante una distribución binomial negativa con parámetro de agregación k (consulte la pág. 315). El principal problema con los recuentos basados en cuadrantes es que dependen en gran medida de la escala. El mismo patrón espacial podría parecer regular cuando se analiza con cuadrantes pequeños, agregado cuando se analiza con cuadrantes de tamaño mediano, pero aleatorio cuando se analiza con cuadrantes grandes.
Las medidas de distancia son de dos tipos generales: medidas desde individuos hasta sus vecinos más cercanos y medidas desde puntos aleatorios hasta el individuo más cercano. Recuerde que el individuo más cercano a un punto aleatorio no es un individuo seleccionado al azar: este protocolo favorece la selección de individuos aislados e individuos en los bordes de los grupos.
En otras circunstancias, es posible que esté dispuesto a dar por sentada la existencia de parches y a realizar un análisis más sofisticado de los atributos espaciales de los parches en sí, su tamaño medio y la variación de tamaño, la variación espacial en el espaciado de parches. de diferentes tamaños, etc.
26.1.1 Puntos aleatorios en un círculo
El círculo está especificado por las coordenadas \(x\) y \(y\) de su centro y por el radio. Podemos calcular las coordenadas de la circunferencia de un círculo de radio r, con su centro ubicado en \(x = y = 0\) así:
#x_c <- r*sin(angle)
#y_c <- r*cos(angle)donde el ángulo varía entre 0 y 2π radianes. Hay dos formas de generar puntos aleatorios dentro de este círculo: una es asumir que el círculo es un objetivo y que estoy apuntando al centro; y el otro es asumir que estoy cortando un parche circular de un mar de puntos espacialmente independientes. En el primer caso, podríamos generar un ángulo aleatorio uniforme, luego generar una distancia aleatoria distribuida uniformemente a lo largo de este radio. En general, estos puntos se agruparán alrededor del centro del círculo porque los radios aleatorios están más densamente agrupados aquí.
point <- function(r) {
angle <- runif(1)*2*pi
length <- runif(1)*r
x_1 <- length*sin(angle)
y_1 <- length*cos(angle)
return (data.frame(x_1,y_1))
}El este y el norte del centro del círculo son e0 y n0 respectivamente, el radio es r y queremos trazar 1000 puntos aleatorios dentro del círculo:
e0 <- 10
n0 <- 10
plot(e0,n0,ylab="",xlab="",ylim=c(0,2*n0),xlim=c(0,2*e0),type="n")
n <- 1000
r <- 10
for (i in 1:n) {
a <- point(r)
e <- e0+a[1]
n <- n0+a[2]
points(e,n,pch=16,col="cyan")
}
Si, en lugar de trazar los puntos aleatorios, dibuja líneas desde el centro del círculo hasta los puntos aleatorios, puede ver exactamente por qué este algoritmo da puntos aleatorios que se agrupan alrededor del centro.
e0 <- 10
n0 <- 10
plot(e0,n0,ylab="",xlab="",ylim=c(0,2*n0),xlim=c(0,2*e0),type="n")
n <- 1000
r <- 10
for (i in 1:n) {
a <- point(r)
e <- e0+a[1]
n <- n0+a[2]
lines(c(e0,e),c(n0,n),col= "red")
}
Una pregunta diferente es el caso del “cortador de galletas”. Si lanzo un cuadrante circular en un mapa espacialmente uniforme de puntos aleatorios, ¿cómo se ve la distribución de mis puntos seleccionados al azar? Aquí está el pseudocódigo:
Hacer un mapa cuadrado de n puntos aleatorios (este y norte uniformes).
Haga un polígono para describir la circunferencia de su cuadrante de muestreo circular.
Coloque el cuadrante en el mapa cuadrado y use la función de maptools para preguntar si todos los puntos del mapa están o no dentro de su círculo (la función se llama point.in.polygon y devuelve un 1 para VERDADERO y un cero para FALSO).
Utilice el vector de salida llamado quería para seleccionar los puntos que están en su círculo.
Aquí está el código R para 10000 puntos aleatorios en una región cuadrada cuyo lado tiene una longitud de 10:
n <- 10000
side <- 10
library(maptools)## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()space <- cbind((runif(n)*side),(runif(n)*side))
plot(space)
circle <- function(e,n,r) {
angle <- seq(0,2*pi,2*pi/360)
x - r*sin(angle)
y - r*cos(angle)
return (cbind((x+e),(y+n)))
}Seleccione los puntos aleatorios en un circulo de radio 1 y centrado en (8,8)
n <- 10000
side <- 16
library(maptools)
space <- cbind((runif(n)*side),(runif(n)*side))
plot(space)
circle <- function(e,n,r) {
angle <- seq(0,2*pi,2*pi/360)
x <- r*sin(angle)
y <- r*cos(angle)
return (cbind((x+e),(y+n)))
}
xc <- 8
yc <- 8
rc <- 1
outline <- circle(xc,yc,rc)
wanted <- point.in.polygon(space[,1],
space[,2],
outline[,1],
outline[,2])
points(space[,1][wanted==1],
space[,2][wanted==1],
col="blue",pch=16)
xc <- 5
yc <- 5
rc <- 2
outline<-circle(xc,yc,rc)
wanted<-point.in.polygon(space[,1],
space[,2],
outline[,1],
outline[,2])
points(space[,1][wanted==1],
space[,2][wanted==1],
col="red",pch=16)
Como se pretendía, no hay agrupación de estos puntos alrededor de los centros de los círculos. Si el círculo representa una pequeña fracción del área total del cuadrado, entonces este método es muy ineficiente.
26.2 - Vecinos cercanos
Supongamos que se nos ha planteado el problema de dibujar líneas para unir los pares de vecinos más cercanos de cualquier conjunto de puntos (\(x, y\)) que se mapean en dos dimensiones. Hay tres pasos para la informática:
necesitamos
calcular la distancia a cada vecino;
identificar la distancia vecina más pequeña para cada individuo;
Utilice estas distancias mínimas para identificar a todos los vecinos más cercanos.
Comenzamos generando una distribución espacial aleatoria de 100 individuos simulando sus coordenadas \(x\) y \(y\) a partir de una distribución de probabilidad uniforme:
x_v <- runif(100)
y_v <- runif(100)
dat_vc <- data.frame(x_v,y_v)El parámetro de gráficos pty = “s” hace que el área de trazado sea cuadrada, como si quisiéramos para este un mapa como este:
ggplot(data=dat_vc,aes(x_v,y_v))+
geom_point(color="brown")+
xlab("x")+
ylab("y")
Calcular las distancias es sencillo: para cada individuo usamos Pitágoras para calcular la distancia a todas las demás plantas. La distancia entre dos puntos con coordenadas (x1, y1) y (x2, y2) es d:
\[d=\sqrt{(y_2-y_1)^2+(x_2-x_1)^2} \]
Se escribe una función para esto:
distance <- function(x1,y1,x2,y2) sqrt((x2-x1)^2+(y2-y1)^2)Ahora recorremos cada individuo \(i\) y calculamos un vector de distancias, d, de todos los demás individuos. La distancia del vecino más cercano es el valor mínimo de \(d\), y la identidad del vecino más cercano, nn, se encuentra usando la función which, which (\(d == min (d [-i])\)), que da el subíndice del mínimo valor de \(d\) (el [-i] es necesario para excluir la distancia 0 que resulta de la distancia del i-ésimo individuo de sí mismo). Aquí está el código completo para calcular las distancias vecinas más cercanas, r, y las identidades, nn, para las 100 individios en el mapa:
r <- numeric(100)
nn <- numeric(100)
d <- numeric(100)
for (i in 1:100) {
for (k in 1:100) d[k] <- distance(x_v[i],y_v[i],x_v[k],y_v[k])
r[i] <- min(d[-i])
nn[i] <- which(d==min(d[-i]))
}Ahora dibuje las distancias dentro del gráfico
par(pty="s")
plot(x_v,y_v,pch=21,bg="red")
for (i in 1:100) lines(c(x_v[i],x_v[nn[i]]),c(y_v[i],y_v[nn[i]]),col="green")
Tenga en cuenta que cuando dos puntos están muy juntos y cada uno es el vecino más cercano del otro, puede parecer que un solo punto no está unido a ningún vecino. La siguiente tarea es averiguar cuántos de los individuos están más cerca del borde del área que de su vecino más cercano. Debido a que los márgenes inferior e izquierdo están en y = 0 y x = 0 respectivamente, la coordenada y de cualquier punto da la distancia desde el borde inferior del área y la coordenada x da la distancia desde el margen izquierdo. Solo necesitamos calcular la distancia de cada individuo desde los márgenes superior y derecho del área:
topd <- 1-y_v
rightd <- 1-x_vAhora usamos la función mínima paralela \(pmin\) para calcular la distancia al margen más cercano para cada uno de los 100 individuos:
edge <- pmin(x_v,y_v,topd,rightd)Finalmente, contamos el número de casos en los que la distancia al borde es menor que la distancia al vecino más cercano:
sum(edge<r)## [1] 18Se identifican estos puntos dentro del mapa, encerrandolos en rojo
plot(x_v,y_v,pch=16)
id <- which(edge<r)
points(x_v[id],y_v[id],col="red",cex=1.5,lwd=2)
Es la distancia vertical u horizontal al borde la que se ha utilizado para identificar estos puntos, por lo que algunos de ellos miran sospechosamente cerca de sus vecinos (por ejemplo, en la esquina inferior izquierda). Los efectos de borde son potencialmente muy importantes en los procesos de puntos espaciales, especialmente cuando hay pocos individuos o el área mapeada es larga y delgada (en lugar de cuadrada o circular). Excluir a los individuos que están más cerca del borde que de su vecino más cercano reduce la distancia media del vecino más cercano:
id <- which(edge<r)
mean(r)## [1] 0.0492494mean(r[-id])## [1] 0.0449517926.2.1 Teselación
El procedimiento de dividir una superficie bidimensional en un mosaico al reducir a la mitad la distancia entre pares de puntos vecinos se llama teselación. Hay una función para hacer esto en el paquete tripack de Albrecht Gebhardt:
x_ts <-runif(100)
y_ts <-runif(100)
dat_ts <- data.frame(x_ts,y_ts) Cree un objeto Voronoi (aquí llamado mapa) aplicando la función llamada voronoi.mosaic a los vectores de las coordenadas \(x\) y \(y\).
mapa <- voronoi.mosaic(x_ts,y_ts)empiece produciendo un diagrama de puntos aleatorios en verde:
ggplot(dat_ts,aes(x_ts,y_ts))+
geom_point(color="green")+
xlab("x")+
ylab("y")
plot(x_ts,y_ts,pch=16,col="green")
plot.voronoi(mapa,pch=16,add=T)
Como puede ver, es relativamente inusual que los puntos estén en el “centro de gravedad” de su parche en mosaico. Cada nodo (círculo negro) es un centro circuncírculo de algún triángulo de la triangulación de Delaunay.
26.3 - Prueba para aleatorización espacial
Clark y Evans (1954) dan una prueba muy simple de aleatoriedad espacial. Suponiendo firmemente que conoce la densidad de población de los individuos, \(\rho\) (generalmente no lo sabe, y necesitaría estimarlo de forma independiente), entonces la distancia media esperada al vecino más cercano es:
\[E(r)=\frac{\sqrt{\rho}}{2}\]
En nuestro ejemplo tenemos 100 individuos en un cuadrado unitario, entonces \(\rho = 0.01\) y \(E (r) = 0.05\). La distancia media real al vecino más cercano fue:
mean(r)## [1] 0.0492494lo cual está muy cerca de la expectativa: esta es claramente una distribución aleatoria de individuos (como la construimos). Un índice de aleatoriedad viene dado por la razón \(\frac{\bar{r}}{E(r)}=\frac{2\bar{r}}{\sqrt{\rho}}\). Esto toma el valor 1 para patrones aleatorios, más de 1 para patrones regulares (espaciados) y menos de 1 para patrones agregados.
Un problema con tales estimaciones de patrón espacial de primer orden (incluidas medidas como la relación varianza-media) es que no pueden dar una idea de la forma en que cambia la distribución espacial dentro de un área.
26.3.1
Las propiedades de segundo orden de un proceso de puntos espaciales describen la forma en que las interacciones espaciales cambian a través del espacio. Estas son medidas computacionalmente intensivas que toman un rango de distancias dentro del área, calculan una medida de patrón, luego trazan un gráfico de la función contra la distancia, para mostrar cómo la medida del patrón cambia con la escala. La medida de segundo orden más utilizada es la función K, que se define como:
\[K(d)=\frac{1}{\lambda}E[número~de~puntos\leq distancia~d~de~un~punto~arbitrario]\]
donde \(\lambda\) es el número medio de puntos por unidad de área (la intensidad del patrón). Si no hay dependencia espacial, entonces el número esperado de puntos que están dentro de una distancia \(d\) de un punto arbitrario es \(\pi d2\) veces la densidad media. Entonces, si la densidad media es de 2 puntos por metro cuadrado (\(\lambda = 2\)), entonces el número esperado de puntos dentro de un radio de 5 m es \(\lambda \pi d2 = 2 × \pi × 52 = 50\pi = 157.1\). Si hay agrupamiento, entonces esperamos un exceso de puntos a distancias cortas (es decir, \(K (d)> \pi d2\) para d pequeña). Del mismo modo, para un patrón espaciado regularmente, esperamos un exceso de distancias largas y, por lo tanto, pocos individuos a distancias cortas (es decir, \(K (d) <\pi d2\)). La K de Ripley (publicada en 1976) se calcula de la siguiente manera:
\[\hat{K}(d)=\frac{1}{n^2}|A|\sum \sum_{i\neq j}{\frac{I_d(d_{ij})}{w_{ij}}}\]
Aquí n es el número de puntos en la región A con área |A|, y \(d_{ij}\) son las distancias entre puntos (la distancia entre los puntos \(i\) y \(j\), para ser precisos). Para tener en cuenta los efectos de borde, el modelo incluye el término \(w_{ij}\), que es la fracción del área, centrada en \(i\) y que pasa por \(j\), que se encuentra dentro del área A (todas las \(w_{ij}\) son 1 para puntos que se encuentran muy lejos de los bordes de la zona). Id (\(d_{ij}\)) es una función indicadora para mostrar qué puntos se deben contar como vecinos en este valor de \(d\): toma el valor 1 si \(d_{ij}\leq d\) y cero en caso contrario (es decir, los puntos con \(d_{ij} > d\) se omiten de la suma). La medida del patrón se obtiene trazando \(\hat{K} (d)\) contra \(d\). Luego, esto se compara con la curva que se observaría bajo una aleatoriedad espacial completa (es decir, una gráfica de \(\pi d2\) contra \(d\)). Cuando se produce la agrupación, \(K (d)> \pi d2\) y la curva se encuentra por encima de la curva de CSR, mientras que los patrones regulares producen una curva por debajo de la curva de CSR. Puede ver por qué necesita la corrección de bordes con este sencillo experimento de simulación. Para el individuo número 1, con coordenadas (x1, y1), calcule las distancias a todos los demás individuos, usando la función \(distance\) que escribimos anteriormente (p. 830):
distances <- numeric(100)
for(i in 1:100) distances[i] <- distance(x_v[1],y_v[1],x_v[i],y_v[i])Ahora averigüe cuántos otros individuos se encuentran a una distancia d de este individuo. Tomemos como ejemplo \(d = 0.1\).
sum(distances<0.1)-1## [1] 6Había otros cuatro individuos dentro de una distancia \(d = 0.1\) del primer individuo (la distancia 0 de sí mismo se incluye en la suma, por lo que tenemos que corregir esto restando 1). El siguiente paso es generalizar el procedimiento de este individuo a todos los individuos. Hacemos una matriz bidimensional llamada dd para contener todas las distancias de cada individuo (filas) a cada otro individuo (columnas):
dd <- numeric(10000)
dd <- matrix(dd,nrow=100)La matriz de distancias se calcula dentro de los bucles tanto para el individuo (\(j\)) como para el vecino (\(i\)) así:
for (j in 1:100) {for(i in 1:100) dd[j,i] <- distance(x_v[j],y_v[j],x_v[i],y_v[i])}Alternativamente, puede usar \(sapply\) con una función anónima como esta, que tiene la ventaja de que no necesitamos preparar la matriz dd por adelantado:
dd <- sapply(1:100,function (i,j=1:100) distance(x_v[j],y_v[j],x_v[i],y_v[i]))Debemos comprobar que el número de individuos dentro de 0.1 del individuo 1 sigue siendo 4 bajo esta nueva notación. Tenga en cuenta el uso de subíndices en blanco [1,] para significar “todas las personas en la fila número 1”:
sum(dd[1,]<0.1)-1## [1] 6Está bien. Queremos calcular la suma de esta cantidad sobre todos los individuos, no solo el número individual
sum(dd<0.1)-100## [1] 336Esto significa que hay 270 casos en los que otros individuos se cuentan dentro de \(d = 0.1\) de los individuos focales. A continuación, cree un vector que contenga un rango de distancias diferentes, \(d\), sobre el cual queremos calcular \(K (d)\) contando el número de individuos dentro de la distancia d, sumados a todos los individuos:
d <- seq(0.01,1,0.01)Para cada una de estas distancias, necesitamos calcular el número total de vecinos de todos los individuos. Entonces, en lugar de 0.1 (en la suma, arriba), necesitamos poner cada uno de los valores \(d\) por turno. El conteo de individuos será un vector de longitud 100 (uno por cada \(d\)):
count <- numeric(100)Calcula la cuenta para cada distancia \(d\):
for (i in 1:100) count[i] <- sum(dd<d[i])-100El recuento esperado aumenta con d como \(\pi d2\), por lo que escalamos nuestro recuento dividiendo por el cuadrado del número total de individuos \(n2 = 1002 = 10000\):
K <- count/10000Finalmente, trace el gráfico \(K\) contra \(d\):
No es sorprendente que cuando muestreamos el área completa (d = 1), contamos a todos los individuos en cada vecindario (K = 1). Para CSR, el gráfico debe seguir πd2, por lo que agregamos una línea para mostrar esto:
plot(d,K,type="l",col="black",ylim=c(0,0.8))
lines(d,pi*d^2,col="blue")
Hasta aproximadamente \(d = 0.2\), la concordancia entre las dos líneas es razonablemente buena, pero para distancias más largas nuestro algoritmo cuenta muy pocos vecinos. Esto se debe a que gran parte del área escaneada alrededor de los individuos marginales es invisible, ya que se encuentra fuera del área de estudio (puede que haya individuos por ahí, pero nunca lo sabremos). Este modelo simple demuestra que la corrección de los bordes es una parte fundamental de la \(K\) de Ripley.
Afortunadamente, no tenemos que escribir una función para calcular un valor corregido para \(K\); está disponible como \(Kfn\) en la biblioteca \(spatial\) incorporada. Aquí lo usamos para analizar el patrón de árboles en el marco de datos llamado \(pines\). La función de biblioteca ppinit lee los datos de un archivo de biblioteca llamado \(pines.dat\) que se almacena en el directorio espacial / ppdata. Luego lo convierte en una lista con los nombres \(\$x\), \(\$y\) y \(\$area\). La primera fila del archivo contiene el número de árboles (71), la segunda fila tiene el nombre del conjunto de datos (\(pines\)), La primera fila del archivo contiene el número de árboles (71), la segunda fila tiene el nombre del conjunto de datos (pinos), la tercera fila tiene los cuatro límites de la región más el factor de escala (0, 96, 0, 100 , 10 de modo que las coordenadas de x inferior y superior se calculan como 0 y 9,6, y las coordenadas de y inferior y superior son 0 y 10). Las filas restantes del archivo de datos contienen las coordenadas xey para cada árbol individual, y estas se convierten en una lista de valores \(x\) y una lista separada de valores \(y\). Necesita conocer estos detalles sobre la estructura de los archivos de datos para poder utilizar estas funciones de biblioteca con sus propios datos (consulte la p. 845).
pines <- ppinit("pines.dat")Primero, configure el área de trazado con dos marcos cuadrados:
windows(7,4)
par(mfrow=c(1,2),pty="s")A la izquierda, haz un mapa usando las ubicaciones xey de los árboles, y a la derecha haz una gráfica de Let) (la medida del patrón) contra la distancia:
par(mfrow=c(1,2),pty="s")
plot(pines,pch=16, col="blue")
plot(Kfn(pines,5),type="s",xlab="distance",ylab="L(t)")
lims <- Kenvl(5,100,Psim(71))
lines(lims$x,lims$lower,lty=2,col="red")
lines(lims$x,lims$upper,lty=2,col="red")
Existe una sugerencia de que a distancias relativamente pequeñas (alrededor de 1 más o menos), los árboles se distribuyen con bastante regularidad (más espaciados que aleatorios), porque la gráfica de L (t) contra la distancia cae por debajo de la envolvente inferior de la línea CSR ( debería estar entre los dos límites en toda su extensión si hubiera RSE). El mecanismo subyacente a esta regularidad espacial (por ejemplo, reclutamiento o mortalidad no aleatorios, competencia entre árboles en crecimiento o ausencia de azar subyacente en el sustrato) debería investigarse en detalle. Con un patrón agregado, la línea caería por encima del sobre superior
Métodos basados en cuadrantes
Otro enfoque para probar la aleatoriedad espacial es contar el número de individuos en cuadrantes de diferentes tamaños. Aquí, los cuadrantes tienen un área de 0.01, por lo que el número esperado por cuadrante es 1. Anteriormente, generamos 100 coordenadas aleatorias para \(x\) y \(y\):
x_q <- runif(100)
y_q <- runif(100)
plot(x_q,y_q,pch=16,col="red")
grid(10,10,lty=1)
Tenga en cuenta que la función \(grid\) no ha hecho exactamente lo que pretendíamos (las cuadrículas no están exactamente en las marcas de verificación). Para contar el número de individuos en cada una de las celdas del mapa, el truco consiste en usar cut para convertir las coordenadas \(x\) y \(y\) del mapa en números de bin (entre 1 y 10 para el tamaño del cuadrante que hemos dibujado aquí). Para lograr esto, los puntos de ruptura son generados por la secuencia (0,1,0.1):
x_t <- cut(x_q,seq(0,1,0.1))
y_t <- cut(y_q,seq(0,1,0.1))Esto crea vectores de subíndices enteros entre 1 y 10 para \(x_t\) y \(y_t\). Ahora todo lo que tenemos que hacer es usar la tabla para contar el número de individuos en cada celda (es decir, en cada combinación de \(x_t\) e \(y_t\)):
count_q <- as.vector (table(x_t,y_t))
table(count_q)## count_q
## 0 1 2 3 4 5 6
## 35 44 13 5 1 1 1Esto muestra que 39 celdas están vacías, 1 celda tenía cinco individuos, pero ninguna celda contenía seis o más individuos. Ahora necesitamos ver cómo se vería esta distribución bajo una hipótesis nula particular.
\[P(x)=\frac{e^{-\lambda}\lambda^x}{x!}\]
Tenga en cuenta que la media depende del tamaño del cuadrante que hemos elegido. Con 100 individuos en toda el área, el número esperado en cualquiera de nuestras 100 celdas, λ, es 1.0. Por lo tanto, las frecuencias esperadas de conteos entre 0 y 5 están dadas por
(expected <- 100*exp(-1)/sapply(0:5,factorial))## [1] 36.7879441 36.7879441 18.3939721 6.1313240 1.5328310 0.3065662El ajuste entre observado y esperado es casi perfecto (como deberíamos esperar, por supuesto, habiendo generado el patrón aleatorio nosotros mismos). En la p. Se muestra una prueba de la significancia de la diferencia entre una distribución de frecuencia observada y esperada.
Patrones agregados y datos de recuento de cuadrantes
A continuación se muestra un ejemplo de un análisis basado en cuadrantes de un patrón espacial agregado. Comenzamos produciendo un mapa de los árboles, luego usamos \(abline\) en lugar \(grid\)) para asegurarnos de que las líneas estén exactamente donde queremos que estén:
trees <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/trees.txt",header = T)
attach(trees)
names(trees)## [1] "x" "y"x_trees <- trees$x
y_trees <- trees$y
plot(x_trees,y_trees,pch=16,col="blue")
abline(v=seq(0,1,0.1),col="lightgray",lty=1)
abline(h=seq(0,1,0.1),col="lightgray",lty=1)
Cortamos los datos y tabulamos los recuentos:
x_t_1 <- cut(x_trees,10)
y_t_1 <- cut(y_trees,10)
count_t_1 <- as.vector(table(x_t_1,y_t_1))
table(count_t_1)## count_t_1
## 0 1 2 3 4 5 6 7 8
## 21 16 15 16 14 10 2 3 3Hay cuadrantes con hasta 8 individuos y, a pesar de que el número medio es mayor a 3 individuos por cuadrado, todavía hay 21 cuadrados completamente vacíos. Las frecuencias esperadas bajo la hipótesis nula de un patrón aleatorio dependen solo del número medio por celda,
mean(count_t_1)## [1] 2.57y como una estimación preliminar de la desviación de la aleatoriedad, calculamos la razón varianza-media (recuerde que con la distribución de Poisson la varianza es igual a la media):
var(count_t_1)/mean(count_t_1)## [1] 1.747082Estos datos están claramente agregados (la relación varianza-media es mayor que 1), por lo que podríamos comparar los conteos con una distribución binomial negativa. Las frecuencias esperadas se estiman multiplicando nuestro número total de cuadrados (100) por las densidades de probabilidad de una distribución binomial negativa generada por la función dnbinom. Esto tiene tres argumentos: los conteos para los que queremos las probabilidades (0:10), la media (\(\mu = 2.57\)) y el parámetro de agregación \(k = \frac{\mu^2}{(var-\mu)} = size = 3.44\):
mean(count_t_1)^2/(var(count_t_1)-mean(count_t_1))## [1] 3.440052Las frecuencias esperadas son:
(expected_t <- sort(dnbinom(0:10, size=3.44, mu=2.57)*100))## [1] 0.3895084 0.7322136 1.3470869 2.4139398 4.1859480 6.9589543
## [7] 10.9366019 14.6700927 15.8853897 20.4862080 21.5799261Estos están razonablemente cerca de las frecuencias observadas (arriba) pero necesitamos cuantificar la falta de ajuste. El plan es mostrar las frecuencias observadas y esperadas como pares de barras, una al lado de la otra
expected_t1 <- expected_t[11:3]
count_t_2 <- c(21,16,15,16,14,10,2,3,3)
dat_exp <- data.frame(count_t_2,expected_t1)
fig_exp <- plot_ly(data=dat_exp,y=~count_t_2,type='bar',name = "Observado")
fig_exp <- fig_exp %>% add_trace(y=~expected_t1,name = "Esperado")
fig_exp <- fig_exp %>% layout(xaxis=list(title=""),
yaxis=list(title="Frecuencias"))
fig_expEl ajuste es razonablemente bueno, pero necesitamos una estimación cuantitativa de la falta de acuerdo entre las distribuciones observadas y esperadas. El chi-cuadrado de Pearson es quizás el más simple. Necesitamos recortar los vectores observados y esperados para que ninguna de las frecuencias esperadas sea menor que 2. La inspección muestra que la frecuencia esperada más baja mayor que 2 está en la ubicación 7, por lo que acumularemos todas las frecuencias en las ubicaciones 6 y superiores
expected_t1[7] <- sum(expected_t1[7:length(expected_t1)])
expected_t1 <- expected_t1[-c(8:length(expected_t1))]
count_t_2[7] <- sum(count_t_2[7:length(count_t_2)])
count_t_2 <- count_t_2[-c(8:length(count_t_2))]Ahora calcule el chi cuadrado de Pearson como \(\sum[(O-E)^2/2]\)
sum((count_t_2-expected_t1)^2/expected_t1)## [1] 3.355271El número de grados de libertad es el número de comparaciones legítimas (7) menos el número de parámetros estimados a partir de los datos (2) menos 1 para la contingencia (es decir, \(7-2-1 = 4\) d.f.). Entonces, la probabilidad de obtener un valor de chi cuadrado de este tamaño (\(3.35\)) o mayor es:
1-pchisq(3.35,4)## [1] 0.5010494Concluimos que el binomio negativo es una descripción imperfecta de estos datos cuadráticos (porque \(p <0.05\)). La razón de la importante falta de ajuste es la grave subestimación de los cuadrantes que contienen solo un árbol y el exceso de cuadrantes que contienen seis o siete árboles.
CVontando cosas en mapas
La convención es que si un punto cae exactamente en el eje x o exactamente en el eje y, entonces se cuenta como si estuviera dentro del área (puntos verdes en el mapa de abajo), pero si cae en el eje superior o en la derecha. eje de la mano, entonces está fuera del área (puntos rojos).
plot(c(0,2),c(0,2),type="n",xlab="",ylab="")
lines(c(0,1,1,0,0),c(0,0,1,1,0))
points(c(0.5,0),c(0,0.5),pch=16,col="green",cex=1.5)
points(c(0.5,1),c(1,0.5),pch=16,col="red",cex=1.5)
De esa manera, todos los puntos en el mapa tienen la misma probabilidad de ser contados y no hay doble conteo de puntos. Los puntos rojos en el eje superior se contarán en el siguiente cuadrante al norte, y los puntos rojos en el eje de la derecha se contarán en el siguiente cuadrante al este.
Esta convención de recuento está incorporada en la función R llamada \(cut\), utilizando corchetes y corchetes como este: “(b1, b2]”, “(b2, b3]”… Para el valor predeterminado derecho = VERDADERO; aquí, el círculo corchete significa ‘mayor que’ y el corchete significa ‘menor o igual que’. Para nuestra convención, necesitamos especificar right = FALSE. Ahora tenemos “[b1, b2)”, “[b2, b3)”. . . donde el corchete significa “mayor o igual que” y el corchete significa “menor que”. A continuación, se muestran algunos datos de prueba:
x_cp <- runif(52,0,1)
length(x_cp)## [1] 52y_cp <- runif(52,0,1)
plot(x_cp,y_cp,pch=16,col="green")
lines(c(0,1,1,0,0),c(0,0,1,1,0),lty=2)
lines(c(0.2,0.2),c(0,1),lty=2)
lines(c(0.4,0.4),c(0,1),lty=2)
lines(c(0.6,0.6),c(0,1),lty=2)
lines(c(0.8,0.8),c(0,1),lty=2)
lines(c(0,1),c(0.2,0.2),lty=2)
lines(c(0,1),c(0.4,0.4),lty=2)
lines(c(0,1),c(0.6,0.6),lty=2)
lines(c(0,1),c(0.8,0.8),lty=2)
Queremos contar el número de puntos en cada uno de los 25 cuadrantes que miden 0,2 × 0,2. La función para lograr esto se corta. Necesitamos especificar seis valores en los que cortar el eje \(x\) y seis valores en los que cortar el eje \(y\). La función convierte el vector de valores continuos de \(x_cp\) en un factor \(x_c\) con cinco niveles (y convierte las coordenadas y en niveles de factor dentro de \(y_c\) de la misma manera):
xc <- cut(x_cp,seq(0,1,0.2),right = FALSE)
yc <- cut(y_cp,seq(0,1,0.2),right = FALSE)Contamos el número de puntos en cada cuadrante usando una tabla, así:
table(yc,xc)## xc
## yc [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1)
## [0,0.2) 2 3 3 2 3
## [0.2,0.4) 3 1 0 2 4
## [0.4,0.6) 2 4 4 0 1
## [0.6,0.8) 1 1 1 3 2
## [0.8,1) 3 1 1 5 0Podemos arreglar esto reordenando los niveles de factor de \(y_c\) usando \(rev\) para invertir el orden de las filas:
yc <- factor(yc,rev(levels(yc)))table(yc,xc)## xc
## yc [0,0.2) [0.2,0.4) [0.4,0.6) [0.6,0.8) [0.8,1)
## [0.8,1) 3 1 1 5 0
## [0.6,0.8) 1 1 1 3 2
## [0.4,0.6) 2 4 4 0 1
## [0.2,0.4) 3 1 0 2 4
## [0,0.2) 2 3 3 2 326.4 - Paquetes para estadística espacial
Además de la biblioteca espacial incorporada, hay dos paquetes de contribución sustancial para analizar datos espaciales. La biblioteca \(Spatstat\) es lo que necesita para el análisis estadístico de patrones de puntos espaciales, mientras que la biblioteca \(spdep\) es buena para el análisis espacial de datos de regiones mapeadas.
Con patrones de puntos, las cosas que querrá hacer incluyen:
Creación, manipulación y trazado de patrones de puntos.
Análisis de datos exploratorios.
Simulación de modelos de procesos puntuales.
Ajuste del modelo paramétrico.
Pruebas de hipótesis y diagnósticos.
mientras que con los mapas podrías
calcular estadísticas espaciales básicas como la I de Moran y la C de Geary.
Crear objetos vecinos de clase nb.
Crear objetos de lista de pesos de clase lw.
Calcular las relaciones vecinas a partir de polígonos (contornos de regiones).
Regiones mapeadas en color sobre la base de estadísticas derivadas.
Debe tomarse un tiempo para dominar los diferentes formatos de los objetos de datos utilizados por los dos paquetes. Perderá mucho tiempo si intenta utilizar las funciones de estas bibliotecas con sus propios archivos de datos no reconstruidos. Aquí está el código para instalar y leer sobre \(Spatstat\) y \(spdep\):
26.4.1 - El paquete spatstat
Necesita usar la función ppp para convertir sus datos de coordenadas en un objeto de clase \(ppp\) que representa un conjunto de datos de patrón de puntos en el plano bidimensional. Nuestro siguiente marco de datos contiene información sobre la ubicación y el tamaño de 3359 plantas de hierba cana en un mapa de 30 m × 15 m:
data_rwm <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/ragwortmap2006.txt",header = T)
attach(data_rwm)
names(data_rwm)## [1] "ROW" "COLUMN" "X" "Y" "rosette"
## [6] "regenerating" "stems" "diameter" "xcoord" "ycoord"
## [11] "type"Las plantas se clasifican como pertenecientes a uno de cuatro tipos: los esqueletos son tallos muertos de plantas que florecieron el año anterior, el rebrote son esqueletos que tienen brotes vivos en la base, las plántulas son plantas pequeñas (de unas pocas semanas) y las rosetas son plantas más grandes. (uno o más años) destinados a florecer este año. La función ppp requiere vectores separados para las coordenadas \(x\) y \(y\): estos están en nuestro archivo con los nombres \(xcoord\) y \(ycoord\). Los argumentos tercero y cuarto de \(ppp\) son las coordenadas de límite para \(x\) y \(y\) respectivamente (en este ejemplo c (0,3000) para \(x\) y c (0,1500) para \(y\)). El argumento final de \(ppp\) contiene un vector de lo que se conoce como “marks”: estos son los niveles de factor asociados con cada uno de los puntos (en este caso, el tipo es esqueleto, rebrote, plántula o roseta). Le da un nombre al objeto \(ppp\) (\(ragwort\)) y lo define así:
xcoord <- as.numeric(data_rwm$xcoord)
ycoord <- as.numeric(data_rwm$ycoord)
type_rgw <- as.factor(data_rwm$type)
ragwort <- ppp(xcoord,ycoord,c(0,3000),c(0,1500),marks=type_rgw)Ahora puede usar el objeto llamado \(ragwort\) en una gran cantidad de funciones diferentes para trazar y modelar estadísticos dentro de la biblioteca Spatstat. Por ejemplo, aquí hay mapas de los patrones de puntos para los cuatro tipos de plantas por separado:
plot(split(ragwort),main="",pch=16)
summary(ragwort)## Marked planar point pattern: 3359 points
## Average intensity 0.0007464444 points per square unit
##
## Coordinates are integers
## i.e. rounded to the nearest unit
##
## Multitype:
## frequency proportion intensity
## regrowth 135 0.04019053 3.000000e-05
## rosette 146 0.04346532 3.244444e-05
## seedling 1100 0.32747840 2.444444e-04
## skeleton 1978 0.58886570 4.395556e-04
##
## Window: rectangle = [0, 3000] x [0, 1500] units
## Window area = 4500000 square unitsque calcula la frecuencia y la intensidad de cada marca (“intensidad” es la densidad media de puntos por unidad de área). En este caso, donde las distancias están en centímetros, la intensidad es el número medio de plantas por centímetro cuadrado (la intensidad más alta son los esqueletos, con 0.000 44 cm – 2). La función \(quadratcount\) produce un resumen útil de recuentos:
plot(quadratcount(ragwort),main="")
Este es el valor predeterminado, pero puede especificar el número de cuadrantes en las direcciones \(x\) y \(y\) (por defecto 5 y 5), o proporcionar vectores numéricos que den las coordenadas xey de los límites de los cuadrantes. Si queremos recuentos en cuadrados de 0,5 m:
plot(quadratcount(ragwort,
xbreaks=c(0,500,1000,1500,2000,2500,3000),
ybreaks=c(0,500,1000,1500)),main="")
Hay funciones para producir gráficos de densidad del patrón de puntos:
Z_rgw <- density.ppp(ragwort)
plot(Z_rgw,main="")
La descripción gráfica clásica de los patrones de puntos espaciales es la K de Ripley
K <- Kest(ragwort)## number of data points exceeds 3000 - computing border correction estimate onlyplot(K, main = "K function")
La línea de puntos roja muestra el número esperado de plantas dentro de un radio r de una planta bajo el supuesto de completa aleatoriedad espacial. La curva observada (negra) se encuentra por encima de esta línea, lo que indica una fuerte agregación espacial en todas las escalas espaciales hasta más de 300 cm.
La función de correlación de pares pcf para los datos de hierba cana se ve así:
pc <- pcf(ragwort)
plot(pc, main = "Pair correlation function")
Existe una fuerte correlación entre pares de plantas a escalas pequeñas, pero mucho menos por encima de r = 20 cm. La función \(distmap\) muestra el mapa de distancias alrededor de plantas individuales:
Z_rgw_1 <- distmap(ragwort)
plot(Z_rgw_1,main="")
Puede utilizar \(Spatstat\) para generar una amplia gama de patrones de puntos aleatorios, incluidos puntos aleatorios uniformes independientes, procesos de puntos de Poisson no homogéneos, procesos de inhibición y procesos de puntos de Gibbs utilizando Metropolis-Hastings (consulte \(?Spatstat\) para obtener más detalles). Algunas funciones útiles sobre distancias de punto a punto en \(Spatstat\) incluyen:
\(nndist\) Distancia del vecino más cercano
\(nnwhich\) Encontrar el vecino más cercano
\(pairdist\) Distancias entre todos los pares de puntos
\(crossdist\) Distacias entre puntos en dos patrones
\(exactdt\) Distancia desde cualquier lugar al punto de información más cercano
\(distmap\) Distancia del mapa
\(density.ppp\) densidad del núcleo suavizada
Hay varias estadísticas de resumen para un patrón de puntos de varios tipos con un componente \(\$marks\) que es un factor:
\(Gcross,~Gdot,~Gmulti\) Multiple distribuciones para los vecinos cercanos.
\(Kcross,~Kdot,~Kmulti\) Multiples funciones K.
\(Jcross,~Jdot,~Jmulti\) Multiples funciones J.
\(Alltypes\) Estimaciones de lo anterior para todos los pares \(i\), \(j\).
\(Lest\) Multiples funciones I.
\(Kcross.inhom\) Contraparte no homogénea de \(Kcross\).
\(Kdot.inhom\) Contraparte no homogénea de \(Kdot\)
Los modelos de procesos puntuales se ajustan utilizando la función ppm de la siguiente manera:
model_36 <- ppm(ragwort, ~marks + polynom(x, y, 2), Poisson())
plot(model_36)
Se han producido ocho mapas de este tipo, que muestran medias (cuatro mapas) y errores estándar (cuatro mapas).
summary(model_36)## Point process model
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.ppp(Q = ragwort, trend = ~marks + polynom(x, y, 2), interaction = Poisson())
## Edge correction: "border"
## [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Marked planar point pattern: 3359 points
## Average intensity 0.000746 points per square unit
## Multitype:
## frequency proportion intensity
## regrowth 135 0.0402 3.00e-05
## rosette 146 0.0435 3.24e-05
## seedling 1100 0.3270 2.44e-04
## skeleton 1980 0.5890 4.40e-04
##
## Window: rectangle = [0, 3000] x [0, 1500] units
## Window area = 4500000 square units
##
## Dummy quadrature points:
## 120 x 120 grid of dummy points, plus 4 corner points
## dummy spacing: 25.0 x 12.5 units
##
## Original dummy parameters: =
## Marked planar point pattern: 67693 points
## Average intensity 0.015 points per square unit
## Multitype:
## frequency proportion intensity
## regrowth 17600 0.260 0.00392
## rosette 17600 0.260 0.00391
## seedling 16700 0.246 0.00370
## skeleton 15800 0.233 0.00351
##
## Window: rectangle = [0, 3000] x [0, 1500] units
## Window area = 4500000 square units
## Quadrature weights:
## (counting weights based on 120 x 120 array of rectangular tiles)
## All weights:
## range: [31.2, 312] total: 1.8e+07
## Weights on data points:
## range: [31.2, 156] total: 424000
## Weights on dummy points:
## range: [31.2, 312] total: 17600000
## --------------------------------------------------------------------------------
## FITTED MODEL:
##
## Nonstationary multitype Poisson process
## Possible marks:
## regrowth rosette seedling skeleton
## ---- Intensity: ----
##
## Log intensity: ~marks + (x + y + I(x^2) + I(x * y) + I(y^2))
##
## Fitted trend coefficients:
## (Intercept) marksrosette marksseedling marksskeleton x
## -1.208572e+01 7.833184e-02 2.097791e+00 2.684567e+00 1.538171e-03
## y I(x^2) I(x * y) I(y^2)
## 1.142238e-03 -3.873534e-07 -1.132803e-07 -3.686329e-07
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -1.208572e+01 1.468160e-01 -1.237348e+01 -1.179797e+01 ***
## marksrosette 7.833184e-02 1.194015e-01 -1.556908e-01 3.123545e-01
## marksseedling 2.097791e+00 9.119484e-02 1.919052e+00 2.276529e+00 ***
## marksskeleton 2.684567e+00 8.895487e-02 2.510218e+00 2.858915e+00 ***
## x 1.538171e-03 1.059553e-04 1.330502e-03 1.745840e-03 ***
## y 1.142238e-03 1.957153e-04 7.586430e-04 1.525833e-03 ***
## I(x^2) -3.873534e-07 2.886002e-08 -4.439180e-07 -3.307888e-07 ***
## I(x * y) -1.132803e-07 5.440746e-08 -2.199170e-07 -6.643681e-09 *
## I(y^2) -3.686329e-07 1.058917e-07 -5.761768e-07 -1.610890e-07 ***
## Zval
## (Intercept) -82.3188462
## marksrosette 0.6560372
## marksseedling 23.0033930
## marksskeleton 30.1789749
## x 14.5171671
## y 5.8362228
## I(x^2) -13.4218001
## I(x * y) -2.0820737
## I(y^2) -3.4812257
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## (Intercept) marksrosette marksseedling marksskeleton x
## -1.208572e+01 7.833184e-02 2.097791e+00 2.684567e+00 1.538171e-03
## y I(x^2) I(x * y) I(y^2)
## 1.142238e-03 -3.873534e-07 -1.132803e-07 -3.686329e-07
##
## Fitted exp(theta):
## (Intercept) marksrosette marksseedling marksskeleton x
## 5.639463e-06 1.081481e+00 8.148148e+00 1.465185e+01 1.001539e+00
## y I(x^2) I(x * y) I(y^2)
## 1.001143e+00 9.999996e-01 9.999999e-01 9.999996e-01produce una tabla enorme de resultados, que incluye lo que los autores denominan “detalles sangrientos”.
26.4.2 - El paquete spdep
La clave para utilizar este paquete es comprender las diferencias entre los distintos formatos en los que se pueden almacenar los datos espaciales:
Coordenadas \(x\) y \(y\) (en una matriz de dos columnas, con x en la columna 1 e y en 2).
Listas de regiones que son vecinas a cada región, con números (potencialmente) desiguales de vecinos en diferentes casos (esto se llama archivo de vecino y pertenece a la clase nb).
Marcos de datos que contienen una región, su vecino y el peso estadístico de la asociación entre las dos regiones en cada fila (class data.frame).
Listas que contienen las identidades de los k vecinos más cercanos (clase knn).
Un objeto de lista de pesos adecuado para calcular la I de Moran o la C de Geary (clase lw);
Listas de polígonos, que definen los contornos de regiones en un mapa (clase polylist).
A diferencia de \(Spatstat\) donde las coordenadas \(x\) y \(y\) estaban en vectores separados, \(spdep\) quiere las coordenadas \(x\) y \(y\) en una única matriz de dos columnas. Para los datos de hierba cana (p. 845) necesitamos escribir:
myco <- cbind(xcoord,ycoord)
myco <- matrix(myco,ncol=2)Una lista sin procesar de coordenadas no contiene información sobre vecinos, pero podemos usar la función \(knearneigh\) para convertir una matriz de coordenadas en un objeto de clase knn. Aquí preguntamos por los cuatro vecinos más cercanos de cada planta:
myco.knn <- knearneigh(myco, k=4)## Warning in knearneigh(myco, k = 4): knearneigh: identical points foundEsta lista de objetos tiene la siguiete estructura.
str(myco.knn)## List of 5
## $ nn : int [1:3359, 1:4] 2 1 4 3 7 4 8 7 10 9 ...
## $ np : int 3359
## $ k : num 4
## $ dimension: int 2
## $ x : num [1:3359, 1:2] 27 29 20 20 78 25 89 97 253 259 ...
## - attr(*, "class")= chr "knn"
## - attr(*, "call")= language knearneigh(x = myco, k = 4)\(\$nn\) contiene 3359 listas, cada una de las cuales es un vector de longitud 4, que contiene las identidades de los cuatro puntos que son los vecinos más cercanos de cada uno de los puntos del 1 al 3359.
\(\$np\) (un número entero) es el número de puntos en el patrón.
\(\$k\) es el número de vecinos de cada punto.
\(\$dimension\) es 2.
\(\$x\) es la matriz de coordenadas de cada punto (\(x\) en la primera columna, \(y\) en la segunda).
Antes de que pueda hacer mucho con un objeto knn, normalmente querrá convertirlo en un objeto vecino (nb) usando la función \(knn2nb\) como esta:
myco.nb <- knn2nb(myco.knn)El concepto esencial para usar el paquete spdep es el objeto vecino (con clase nb). Para una ubicación dada, típicamente identificada por las coordenadas (x, y) de su centroide, el objeto vecino es una lista, con los elementos de la lista numerados del 1 al número de ubicaciones, y cada elemento de la lista contiene un vector. de enteros que representan las identidades de las ubicaciones que comparten un límite con esa ubicación. El punto importante es que es probable que los diferentes vectores tengan diferentes longitudes. Puedes hacer cosas interesantes con nb objects. Aquí hay una gráfica con cada punto unido a sus cuatro vecinos más cercanos: usted especifica el objeto nb y la matriz de coordenadas:
plot(myco.nb,myco,pch=16)
La forma más sencilla de crear un objeto nb es leer un archivo de texto que contenga una fila para cada relación de vecino, utilizando la función de entrada especial read.gwt2nb. La fila de encabezado puede adoptar una de dos formas. El más simple (llamado “GWT de estilo antiguo”) es un único entero que indica el número de ubicaciones en el archivo. Siempre habrá muchas más filas en el archivo de datos que este número, porque cada ubicación normalmente tendrá varios vecinos. La segunda forma de la fila de encabezado tiene cuatro elementos: el primero se establece arbitrariamente en cero, el segundo es el número entero de ubicaciones, el tercero es el nombre del objeto de forma y el cuarto es el vector de nombres que identifican las ubicaciones. Un ejemplo debería aclarar esto. Estos son los contenidos de un texto archivo llamado naydf.txt:
El 5 de la primera fila indica que este archivo contiene información sobre cinco ubicaciones. En las líneas siguientes, el primer número identifica la ubicación, el segundo número identifica a uno de sus vecinos y el tercer número es el peso de esa relación. Por lo tanto, la ubicación 5 tiene solo dos vecinos, y son las ubicaciones 3 y 4 (las dos últimas filas del archivo). Creamos un objeto vecino para estos datos con la función read.gwt2nb como esta:
dd <- read.gwt2nb("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/naydf.txt")## Warning in read.gwt2nb("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/
## naydf.txt"): Old-style GWT fileA continuación, se muestra un resumen del objeto vecino recién creado llamado dd:
summary(dd)## Neighbour list object:
## Number of regions: 5
## Number of nonzero links: 13
## Percentage nonzero weights: 52
## Average number of links: 2.6
## Non-symmetric neighbours list
## Link number distribution:
##
## 2 3
## 2 3
## 2 least connected regions:
## 1 5 with 2 links
## 3 most connected regions:
## 2 3 4 with 3 linksAquí hay 5 vectores de vecinos
dd[[1]]## [1] 2 3dd[[2]]## [1] 1 3 4dd[[3]]## [1] 1 2 5dd[[4]]## [1] 2 3 5dd[[5]]## [1] 3 4las coordenadas de las 5 locaiones necesitan ser especificadas
coox <- c(1,2,3,4,5)
cooy <- c(3,1,2,0,3)y los vectores de coordenadas deben combinarse en una matriz de dos columnas. Ahora podemos usar \(plot\) con \(dd\) y la matriz de coordenadas para indicar las relaciones vecinas de las cinco ubicaciones de esta manera:
plot(dd,matrix(cbind(coox,cooy),ncol=2))
text(coox,cooy,as.character(1:5),pos=rep(3,5))
Tenga en cuenta el uso de \(pos = 3\) para colocar los números de ubicación del 1 al 5 sobre sus puntos. Puede ver que las ubicaciones 1 y 5 son las menos conectadas (dos vecinos) y la ubicación 3 es la más conectada (cuatro vecinos). Tenga en cuenta que la especificación en el archivo de datos no fue completamente recíproca, porque la ubicación 4 se definió como un vecino de ubicación 3 pero no al revés. Hay un comentario, Lista de vecinos no simétricos, en la salida del resumen (dd) para llamar la atención sobre esto. Una función \(make.sym.nb (dd)\) está disponible para convertir el objeto dd en una lista de vecinos simétricos.
Para calcular índices como la I de Moran y la C de Geary, necesita un objeto de “lista de pesos”. Esto se crea más simplemente a partir de un objeto vecino usando la función \(nb2listw\). Para los datos de \(ragwort\), ya hemos creado un objeto vecino llamado \(myco.nb\) y creamos el objeto de lista de pesos myco.lw así:
myco.lw <- nb2listw(myco.nb, style="W")
myco.lw## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 3359
## Number of nonzero links: 13436
## Percentage nonzero weights: 0.1190831
## Average number of links: 4
## Non-symmetric neighbours list
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 3359 11282881 3359 1458.5 14076.5Ahí hay 3 pruebas clásicas basadas en productos espaciales cruzados \(C(i,j)\) donde \(z(i)=[x(i)-\frac{\mu(x)}{sd(x)}]\)
Moran \(C(i,j)=z(i)z(j)\)
Geary \(C(i,j)=(z(i)-z(j))^2\)
Sokal \(C(i,j)=|z(i)-z(j)|\)
Aquí está la prueba del índice de Moran para \(ragwort\), usando el objeto de la lista de pesos \(myco.lw\):
moran(1:3359,myco.lw,length(myco.nb),Szero(myco.lw))## $I
## [1] 0.9931203
##
## $K
## [1] 1.8Aquí está la C de Gary para la misma información
geary(1:3359,myco.lw,length(myco.nb),length(myco.nb)-1,Szero(myco.lw))## $C
## [1] 0.004546136
##
## $K
## [1] 1.8La prueba de permutación de Mantel:
sp.mantel.mc(1:3359,myco.lw,nsim=99)##
## Mantel permutation test for moran measure
##
## data: 1:3359
## weights: myco.lw
## number of simulations + 1: 100
##
## statistic = 3334.9, observed rank = 100, p-value = 0.01
## alternative hypothesis: greater
## sample estimates:
## mean of permutations sd of permutations
## -5.529673 38.517995En todos los casos, el primer argumento es un vector de números de ubicación (1 a 3359 en el ejemplo \(ragwort\)), el segundo argumento es el objeto de lista de peso myco.lw. Para moran, el tercer argumento es la longitud del objeto vecino, length (\(myco.nb\)) y el cuarto es Szero (\(myco.lw4\)), la suma global de pesos, los cuales se evalúan a 3359 en este caso. La función geary tiene un argumento adicional, \(length (myco.nb)\) -1, y \(sp.mantel.mc\) especifica el número de simulaciones.
26.4.3 - Lista de poligonos
Quizás los datos espaciales más complejos manejados por \(spdep\) comprenden contornos digitalizados (conjuntos de coordenadas \(x\) y \(y\)) que definen múltiples regiones, cada una de las cuales puede ser interpretada por R como un polígono. Aquí hay una lista de este tipo del conjunto de datos integrado de \(columbus\):
data(columbus)
polys## [[1]]
## [,1] [,2]
## [1,] 8.624129 14.23698
## [2,] 8.559700 14.74245
## [3,] 8.809452 14.73443
## [4,] 8.808413 14.63652
## [5,] 8.919305 14.63850
## [6,] 9.087138 14.63049
## [7,] 9.099965 14.24483
## [8,] 9.015047 14.24184
## [9,] 9.008951 13.99506
## [10,] 8.818140 14.00205
## [11,] 8.653305 14.00809
## [12,] 8.642902 14.08971
## [13,] 8.632592 14.17059
## [14,] 8.625826 14.22367
## [15,] 8.624129 14.23698
## attr(,"bbox")
## [1] 8.559700 13.995060 9.099965 14.742450
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 15
##
##
## [[2]]
## [,1] [,2]
## [1,] 8.252790 14.23694
## [2,] 8.282758 14.22994
## [3,] 8.330711 14.22994
## [4,] 8.383658 14.22893
## [5,] 8.444600 14.22892
## [6,] 8.544504 14.23490
## [7,] 8.624129 14.23698
## [8,] 8.625826 14.22367
## [9,] 8.632592 14.17059
## [10,] 8.642902 14.08971
## [11,] 8.653305 14.00809
## [12,] 8.662189 13.90990
## [13,] 8.666551 13.86170
## [14,] 8.605282 13.83925
## [15,] 8.579310 13.84125
## [16,] 8.562577 13.84254
## [17,] 8.540359 13.84240
## [18,] 8.516387 13.84168
## [19,] 8.502935 13.83873
## [20,] 8.473408 13.83227
## [21,] 8.459499 13.82035
## [22,] 8.431433 13.79331
## [23,] 8.415447 13.79031
## [24,] 8.387156 13.78897
## [25,] 8.373487 13.78832
## [26,] 8.323546 13.78608
## [27,] 8.284572 13.78433
## [28,] 8.291548 13.74137
## [29,] 8.229603 13.72739
## [30,] 8.226613 13.74438
## [31,] 8.215644 13.79433
## [32,] 8.198687 13.85828
## [33,] 8.169725 13.88326
## [34,] 8.127770 13.89226
## [35,] 8.093802 13.89126
## [36,] 8.063838 13.90526
## [37,] 8.044872 13.94322
## [38,] 8.037889 13.96720
## [39,] 7.999116 14.02457
## [40,] 7.999366 14.03497
## [41,] 8.003014 14.18702
## [42,] 7.950089 14.24397
## [43,] 8.111939 14.26393
## [44,] 8.147892 14.23296
## [45,] 8.181855 14.22596
## [46,] 8.209828 14.22695
## [47,] 8.252790 14.23694
## attr(,"bbox")
## [1] 7.950089 13.727390 8.666551 14.263930
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 47
##
##
## [[3]]
## [,1] [,2]
## [1,] 8.653305 14.00809
## [2,] 8.818140 14.00205
## [3,] 9.008951 13.99506
## [4,] 9.008928 13.93811
## [5,] 9.343591 13.91308
## [6,] 9.351485 13.67529
## [7,] 9.298501 13.58938
## [8,] 9.273822 13.58856
## [9,] 9.244555 13.59138
## [10,] 9.242543 13.55841
## [11,] 9.196582 13.54443
## [12,] 9.190605 13.58639
## [13,] 9.166627 13.58140
## [14,] 9.161684 13.70829
## [15,] 8.909940 13.71553
## [16,] 8.677577 13.72221
## [17,] 8.666551 13.86170
## [18,] 8.662189 13.90990
## [19,] 8.653305 14.00809
## attr(,"bbox")
## [1] 8.653305 13.544430 9.351485 14.008090
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 19
##
##
## [[4]]
## [,1] [,2]
## [1,] 8.459499 13.82035
## [2,] 8.473408 13.83227
## [3,] 8.502935 13.83873
## [4,] 8.516387 13.84168
## [5,] 8.540359 13.84240
## [6,] 8.562577 13.84254
## [7,] 8.579310 13.84125
## [8,] 8.605282 13.83925
## [9,] 8.666551 13.86170
## [10,] 8.677577 13.72221
## [11,] 8.685274 13.63952
## [12,] 8.628179 13.63946
## [13,] 8.588302 13.64182
## [14,] 8.571110 13.64127
## [15,] 8.547516 13.64382
## [16,] 8.537268 13.64443
## [17,] 8.505299 13.64443
## [18,] 8.459344 13.64443
## [19,] 8.450570 13.60454
## [20,] 8.439336 13.60552
## [21,] 8.380410 13.61647
## [22,] 8.385412 13.63445
## [23,] 8.316472 13.61648
## [24,] 8.294443 13.60457
## [25,] 8.279500 13.59650
## [26,] 8.247527 13.58651
## [27,] 8.201574 13.59151
## [28,] 8.201584 13.61449
## [29,] 8.198595 13.63548
## [30,] 8.233589 13.70341
## [31,] 8.229603 13.72739
## [32,] 8.291548 13.74137
## [33,] 8.284572 13.78433
## [34,] 8.323546 13.78608
## [35,] 8.373487 13.78832
## [36,] 8.387156 13.78897
## [37,] 8.415447 13.79031
## [38,] 8.431433 13.79331
## [39,] 8.459499 13.82035
## attr(,"bbox")
## [1] 8.198595 13.586510 8.685274 13.861700
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 39
##
##
## [[5]]
## [,1] [,2]
## [1,] 8.685274 13.63952
## [2,] 8.677577 13.72221
## [3,] 8.909940 13.71553
## [4,] 9.161684 13.70829
## [5,] 9.166627 13.58140
## [6,] 9.190605 13.58639
## [7,] 9.196582 13.54443
## [8,] 9.242543 13.55841
## [9,] 9.244555 13.59138
## [10,] 9.273822 13.58856
## [11,] 9.298501 13.58938
## [12,] 9.310472 13.54541
## [13,] 9.401384 13.55040
## [14,] 9.333297 13.27242
## [15,] 9.236267 12.87628
## [16,] 9.233387 12.86400
## [17,] 8.943573 12.86222
## [18,] 8.757729 12.86109
## [19,] 8.733970 13.11634
## [20,] 8.685274 13.63952
## attr(,"bbox")
## [1] 8.677577 12.861090 9.401384 13.722210
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 20
##
##
## [[6]]
## [,1] [,2]
## [1,] 9.401384 13.55040
## [2,] 9.434411 13.69427
## [3,] 9.605247 13.69824
## [4,] 9.651198 13.69224
## [5,] 9.687166 13.69723
## [6,] 9.686146 13.64528
## [7,] 9.845992 13.65225
## [8,] 10.050790 13.65023
## [9,] 10.103720 13.60326
## [10,] 10.175630 13.56529
## [11,] 10.180600 13.48236
## [12,] 10.167600 13.47137
## [13,] 10.153610 13.45439
## [14,] 10.135620 13.43940
## [15,] 10.119630 13.42941
## [16,] 10.121600 13.34449
## [17,] 10.096620 13.34249
## [18,] 10.085630 13.33350
## [19,] 10.052660 13.33650
## [20,] 10.027670 13.29854
## [21,] 9.772106 13.29211
## [22,] 9.677010 13.29659
## [23,] 9.671007 13.27361
## [24,] 9.333297 13.27242
## [25,] 9.401384 13.55040
## attr(,"bbox")
## [1] 9.333297 13.272420 10.180600 13.698240
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 25
##
##
## [[7]]
## [,1] [,2]
## [1,] 8.037741 13.60752
## [2,] 8.062716 13.60452
## [3,] 8.072695 13.58054
## [4,] 8.115653 13.57954
## [5,] 8.130637 13.57654
## [6,] 8.197452 13.57577
## [7,] 8.189590 13.54576
## [8,] 8.179450 13.52983
## [9,] 8.164705 13.51106
## [10,] 8.140126 13.47977
## [11,] 8.128536 13.47010
## [12,] 8.121718 13.46104
## [13,] 8.109848 13.43264
## [14,] 8.104980 13.42099
## [15,] 8.181537 13.41591
## [16,] 8.293398 13.40849
## [17,] 8.380322 13.40466
## [18,] 8.456573 13.10407
## [19,] 8.425150 13.09392
## [20,] 8.412152 13.06895
## [21,] 8.351270 13.06650
## [22,] 8.288272 13.06397
## [23,] 8.280263 13.02600
## [24,] 8.232307 13.02001
## [25,] 8.220303 12.97905
## [26,] 8.154367 12.97806
## [27,] 8.145204 12.94202
## [28,] 8.104401 12.94305
## [29,] 8.062443 12.94410
## [30,] 8.052461 12.96309
## [31,] 8.052497 13.05201
## [32,] 8.048532 13.12994
## [33,] 8.032543 13.11795
## [34,] 8.013560 13.11496
## [35,] 7.989584 13.11596
## [36,] 7.962620 13.13794
## [37,] 7.923680 13.19190
## [38,] 7.898726 13.24386
## [39,] 7.887743 13.25785
## [40,] 7.871759 13.25885
## [41,] 7.868795 13.33878
## [42,] 7.866810 13.37175
## [43,] 7.851837 13.40073
## [44,] 7.844857 13.43270
## [45,] 7.840892 13.50863
## [46,] 7.824926 13.55260
## [47,] 7.803965 13.59656
## [48,] 7.801974 13.61554
## [49,] 7.901874 13.60854
## [50,] 7.902888 13.64451
## [51,] 7.996794 13.63950
## [52,] 7.998781 13.61452
## [53,] 8.037741 13.60752
## attr(,"bbox")
## [1] 7.801974 12.942020 8.456573 13.644510
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 53
##
##
## [[8]]
## [,1] [,2]
## [1,] 8.247527 13.58651
## [2,] 8.279500 13.59650
## [3,] 8.294443 13.60457
## [4,] 8.316472 13.61648
## [5,] 8.385412 13.63445
## [6,] 8.380410 13.61647
## [7,] 8.439336 13.60552
## [8,] 8.450570 13.60454
## [9,] 8.459344 13.64443
## [10,] 8.505299 13.64443
## [11,] 8.537268 13.64443
## [12,] 8.547516 13.64382
## [13,] 8.571110 13.64127
## [14,] 8.588302 13.64182
## [15,] 8.628179 13.63946
## [16,] 8.685274 13.63952
## [17,] 8.733970 13.11634
## [18,] 8.651877 13.11367
## [19,] 8.596990 13.11189
## [20,] 8.585998 13.10689
## [21,] 8.517033 13.10539
## [22,] 8.456573 13.10407
## [23,] 8.380322 13.40466
## [24,] 8.293398 13.40849
## [25,] 8.181537 13.41591
## [26,] 8.104980 13.42099
## [27,] 8.109848 13.43264
## [28,] 8.121718 13.46104
## [29,] 8.128536 13.47010
## [30,] 8.140126 13.47977
## [31,] 8.164705 13.51106
## [32,] 8.179450 13.52983
## [33,] 8.189590 13.54576
## [34,] 8.197452 13.57577
## [35,] 8.201574 13.59151
## [36,] 8.247527 13.58651
## attr(,"bbox")
## [1] 8.10498 13.10407 8.73397 13.64443
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 36
##
##
## [[9]]
## [,1] [,2]
## [1,] 9.333297 13.27242
## [2,] 9.671007 13.27361
## [3,] 9.677010 13.29659
## [4,] 9.772106 13.29211
## [5,] 10.027670 13.29854
## [6,] 9.950614 12.98383
## [7,] 10.004560 12.97683
## [8,] 10.082510 13.03377
## [9,] 10.083650 13.01987
## [10,] 10.086820 12.98143
## [11,] 10.089170 12.95281
## [12,] 10.091630 12.92304
## [13,] 10.093510 12.90022
## [14,] 10.095430 12.87690
## [15,] 10.015440 12.72405
## [16,] 9.763668 12.67313
## [17,] 9.723675 12.59520
## [18,] 9.555828 12.59519
## [19,] 9.471498 12.59571
## [20,] 9.386005 12.59624
## [21,] 9.383021 12.62721
## [22,] 9.258048 12.63061
## [23,] 9.124277 12.63424
## [24,] 9.146788 12.65874
## [25,] 9.166311 12.67998
## [26,] 9.187246 12.70917
## [27,] 9.206208 12.75884
## [28,] 9.213532 12.77803
## [29,] 9.220885 12.80572
## [30,] 9.233387 12.86400
## [31,] 9.236267 12.87628
## [32,] 9.333297 13.27242
## attr(,"bbox")
## [1] 9.124277 12.595190 10.095430 13.298540
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 32
##
##
## [[10]]
## [,1] [,2]
## [1,] 10.08251 13.03377
## [2,] 10.09250 13.05275
## [3,] 10.12649 13.09071
## [4,] 10.18246 13.15765
## [5,] 10.26742 13.25455
## [6,] 10.28939 13.24755
## [7,] 10.31636 13.24455
## [8,] 10.33235 13.23955
## [9,] 10.36826 13.24643
## [10,] 10.39429 13.24654
## [11,] 10.41727 13.24754
## [12,] 10.45024 13.25153
## [13,] 10.47522 13.26152
## [14,] 10.47522 13.27251
## [15,] 10.53416 13.27150
## [16,] 10.56613 13.27150
## [17,] 10.60309 13.26550
## [18,] 10.63206 13.26349
## [19,] 10.63905 13.24451
## [20,] 10.64004 13.21953
## [21,] 10.64902 13.19955
## [22,] 10.64601 13.17657
## [23,] 10.64701 13.15759
## [24,] 10.64600 13.14460
## [25,] 10.63800 13.12961
## [26,] 10.63899 13.10463
## [27,] 10.63198 13.07066
## [28,] 10.62698 13.05168
## [29,] 10.62797 13.02670
## [30,] 10.63196 13.00872
## [31,] 10.63395 12.98374
## [32,] 10.62893 12.94578
## [33,] 10.63991 12.91880
## [34,] 10.63790 12.89082
## [35,] 10.64589 12.86584
## [36,] 10.64787 12.84286
## [37,] 10.64968 12.83018
## [38,] 10.50110 12.80532
## [39,] 10.35660 12.78114
## [40,] 10.17841 12.75132
## [41,] 10.01544 12.72405
## [42,] 10.09543 12.87690
## [43,] 10.09351 12.90022
## [44,] 10.09163 12.92304
## [45,] 10.08917 12.95281
## [46,] 10.08682 12.98143
## [47,] 10.08365 13.01987
## [48,] 10.08251 13.03377
## attr(,"bbox")
## [1] 10.01544 12.72405 10.64968 13.27251
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 48
##
##
## [[11]]
## [,1] [,2]
## [1,] 8.585998 13.10689
## [2,] 8.596990 13.11189
## [3,] 8.651877 13.11367
## [4,] 8.733970 13.11634
## [5,] 8.757729 12.86109
## [6,] 8.700783 12.85710
## [7,] 8.693771 12.81114
## [8,] 8.632830 12.81015
## [9,] 8.629889 12.94503
## [10,] 8.572947 12.95003
## [11,] 8.585998 13.10689
## attr(,"bbox")
## [1] 8.572947 12.810150 8.757729 13.116340
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 11
##
##
## [[12]]
## [,1] [,2]
## [1,] 8.456573 13.10407
## [2,] 8.517033 13.10539
## [3,] 8.585998 13.10689
## [4,] 8.572947 12.95003
## [5,] 8.629889 12.94503
## [6,] 8.632830 12.81015
## [7,] 8.582890 12.80977
## [8,] 8.509166 12.80920
## [9,] 8.487373 12.93066
## [10,] 8.456573 13.10407
## attr(,"bbox")
## [1] 8.456573 12.809200 8.632830 13.106890
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 10
##
##
## [[13]]
## [,1] [,2]
## [1,] 8.145204 12.94202
## [2,] 8.154367 12.97806
## [3,] 8.220303 12.97905
## [4,] 8.232307 13.02001
## [5,] 8.280263 13.02600
## [6,] 8.288272 13.06397
## [7,] 8.351270 13.06650
## [8,] 8.412152 13.06895
## [9,] 8.425150 13.09392
## [10,] 8.456573 13.10407
## [11,] 8.487373 12.93066
## [12,] 8.456042 12.93171
## [13,] 8.402055 12.93353
## [14,] 8.337126 12.93572
## [15,] 8.245318 12.93882
## [16,] 8.145204 12.94202
## attr(,"bbox")
## [1] 8.145204 12.930660 8.487373 13.104070
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 16
##
##
## [[14]]
## [,1] [,2]
## [1,] 8.062443 12.94410
## [2,] 8.104401 12.94305
## [3,] 8.145204 12.94202
## [4,] 8.245318 12.93882
## [5,] 8.337126 12.93572
## [6,] 8.402055 12.93353
## [7,] 8.456042 12.93171
## [8,] 8.487373 12.93066
## [9,] 8.509166 12.80920
## [10,] 8.512938 12.78818
## [11,] 8.478971 12.78819
## [12,] 8.431012 12.78806
## [13,] 8.259181 12.78760
## [14,] 8.122319 12.78723
## [15,] 8.062443 12.94410
## attr(,"bbox")
## [1] 8.062443 12.787230 8.512938 12.944100
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 15
##
##
## [[15]]
## [,1] [,2]
## [1,] 8.757729 12.86109
## [2,] 8.943573 12.86222
## [3,] 9.233387 12.86400
## [4,] 9.220885 12.80572
## [5,] 9.213532 12.77803
## [6,] 9.206208 12.75884
## [7,] 9.187246 12.70917
## [8,] 9.166311 12.67998
## [9,] 9.146788 12.65874
## [10,] 9.124277 12.63424
## [11,] 8.913297 12.60980
## [12,] 8.855528 12.60630
## [13,] 8.785566 12.53237
## [14,] 8.777865 12.62803
## [15,] 8.757729 12.86109
## attr(,"bbox")
## [1] 8.757729 12.532370 9.233387 12.864000
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 15
##
##
## [[16]]
## [,1] [,2]
## [1,] 8.632830 12.81015
## [2,] 8.693771 12.81114
## [3,] 8.700783 12.85710
## [4,] 8.757729 12.86109
## [5,] 8.777865 12.62803
## [6,] 8.785566 12.53237
## [7,] 8.758386 12.49685
## [8,] 8.726524 12.45521
## [9,] 8.709477 12.42731
## [10,] 8.694607 12.41749
## [11,] 8.657633 12.39151
## [12,] 8.605762 12.36794
## [13,] 8.589687 12.36155
## [14,] 8.512938 12.78818
## [15,] 8.509166 12.80920
## [16,] 8.582890 12.80977
## [17,] 8.632830 12.81015
## attr(,"bbox")
## [1] 8.509166 12.361550 8.785566 12.861090
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 17
##
##
## [[17]]
## [,1] [,2]
## [1,] 10.35660 12.78114
## [2,] 10.50110 12.80532
## [3,] 10.64968 12.83018
## [4,] 10.70282 12.83886
## [5,] 10.70979 12.77491
## [6,] 10.66482 12.76293
## [7,] 10.66980 12.70398
## [8,] 10.66777 12.64803
## [9,] 10.66697 12.59503
## [10,] 10.66734 12.58011
## [11,] 10.66970 12.48318
## [12,] 10.42175 12.44026
## [13,] 10.42099 12.57813
## [14,] 10.35606 12.59812
## [15,] 10.35639 12.71201
## [16,] 10.35660 12.78114
## attr(,"bbox")
## [1] 10.35606 12.44026 10.70979 12.83886
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 16
##
##
## [[18]]
## [,1] [,2]
## [1,] 8.431012 12.78806
## [2,] 8.478971 12.78819
## [3,] 8.512938 12.78818
## [4,] 8.589687 12.36155
## [5,] 8.565508 12.35533
## [6,] 8.554715 12.35533
## [7,] 8.523073 12.35531
## [8,] 8.502639 12.35708
## [9,] 8.481381 12.36291
## [10,] 8.438839 12.37356
## [11,] 8.406447 12.39673
## [12,] 8.391569 12.40738
## [13,] 8.369997 12.42818
## [14,] 8.358474 12.44566
## [15,] 8.423884 12.44849
## [16,] 8.431012 12.78806
## attr(,"bbox")
## [1] 8.358474 12.355310 8.589687 12.788190
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 16
##
##
## [[19]]
## [,1] [,2]
## [1,] 8.122319 12.78723
## [2,] 8.259181 12.78760
## [3,] 8.431012 12.78806
## [4,] 8.423884 12.44849
## [5,] 8.358474 12.44566
## [6,] 8.350509 12.48193
## [7,] 8.347977 12.49346
## [8,] 8.348999 12.54742
## [9,] 8.332018 12.54742
## [10,] 8.305043 12.55042
## [11,] 8.301018 12.59309
## [12,] 8.299444 12.60978
## [13,] 8.296090 12.64534
## [14,] 8.223163 12.64534
## [15,] 8.148235 12.64535
## [16,] 8.138268 12.70331
## [17,] 8.138283 12.73728
## [18,] 8.122319 12.78723
## attr(,"bbox")
## [1] 8.122319 12.445660 8.431012 12.788060
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 18
##
##
## [[20]]
## [,1] [,2]
## [1,] 10.015440 12.72405
## [2,] 10.178410 12.75132
## [3,] 10.356600 12.78114
## [4,] 10.356390 12.71201
## [5,] 10.356060 12.59812
## [6,] 10.420990 12.57813
## [7,] 10.421750 12.44026
## [8,] 10.424250 11.99064
## [9,] 10.425640 11.74186
## [10,] 10.252800 11.74282
## [11,] 10.056590 11.74391
## [12,] 10.049910 11.76295
## [13,] 10.049750 11.77896
## [14,] 10.049290 11.82283
## [15,] 10.048910 11.85880
## [16,] 10.042080 11.91276
## [17,] 10.024220 11.96231
## [18,] 10.007160 12.00968
## [19,] 9.991961 12.03151
## [20,] 9.956293 12.08275
## [21,] 9.914780 12.14239
## [22,] 9.899952 12.16369
## [23,] 9.892879 12.17385
## [24,] 9.884351 12.19153
## [25,] 10.001230 12.17054
## [26,] 9.966332 12.34139
## [27,] 9.841084 12.31474
## [28,] 9.935688 12.50187
## [29,] 9.941194 12.51213
## [30,] 9.933469 12.59941
## [31,] 10.015440 12.72405
## attr(,"bbox")
## [1] 9.841084 11.741860 10.425640 12.781140
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 31
##
##
## [[21]]
## [,1] [,2]
## [1,] 7.842435 12.40560
## [2,] 8.131140 12.37360
## [3,] 8.136807 12.34855
## [4,] 8.141801 12.32647
## [5,] 8.148350 12.29872
## [6,] 8.157790 12.28122
## [7,] 8.174397 12.25043
## [8,] 8.190520 12.23764
## [9,] 8.225327 12.20686
## [10,] 8.250281 12.19443
## [11,] 8.310225 12.15693
## [12,] 8.326121 12.15328
## [13,] 8.358864 12.14595
## [14,] 8.408770 12.13477
## [15,] 8.434194 12.13921
## [16,] 8.471713 12.14575
## [17,] 8.495033 12.13777
## [18,] 8.505514 12.12821
## [19,] 8.519205 12.11523
## [20,] 8.538624 12.08879
## [21,] 8.544462 12.07360
## [22,] 8.563573 12.02385
## [23,] 8.525978 11.92466
## [24,] 8.517660 11.91471
## [25,] 8.471678 11.86373
## [26,] 8.455998 11.84761
## [27,] 8.412508 11.80735
## [28,] 8.404299 11.80047
## [29,] 8.386433 11.78549
## [30,] 8.373654 11.77409
## [31,] 8.282732 11.74413
## [32,] 8.256792 11.82806
## [33,] 8.103924 11.78512
## [34,] 8.120859 11.67022
## [35,] 7.733179 11.52739
## [36,] 7.688279 11.66328
## [37,] 7.717283 11.74121
## [38,] 7.580475 11.88210
## [39,] 7.637434 11.91706
## [40,] 7.616512 12.05794
## [41,] 7.403729 12.08095
## [42,] 7.404763 12.16487
## [43,] 7.323862 12.21184
## [44,] 7.196021 12.29578
## [45,] 7.135107 12.35973
## [46,] 7.108153 12.40670
## [47,] 7.092213 12.51261
## [48,] 7.061330 12.72542
## [49,] 7.143250 12.72441
## [50,] 7.155218 12.67446
## [51,] 7.173169 12.59852
## [52,] 7.206114 12.54456
## [53,] 7.270038 12.51058
## [54,] 7.341958 12.48660
## [55,] 7.429863 12.46660
## [56,] 7.459833 12.46360
## [57,] 7.481781 12.38966
## [58,] 7.498770 12.40365
## [59,] 7.516755 12.40964
## [60,] 7.548725 12.41263
## [61,] 7.559715 12.41463
## [62,] 7.558692 12.35668
## [63,] 7.569672 12.33470
## [64,] 7.616611 12.29873
## [65,] 7.669571 12.32769
## [66,] 7.721513 12.30771
## [67,] 7.708548 12.36266
## [68,] 7.639627 12.38964
## [69,] 7.646636 12.43061
## [70,] 7.842435 12.40560
## attr(,"bbox")
## [1] 7.061330 11.527390 8.563573 12.725420
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 70
##
##
## [[22]]
## [,1] [,2]
## [1,] 9.386005 12.59624
## [2,] 9.471498 12.59571
## [3,] 9.555828 12.59519
## [4,] 9.723675 12.59520
## [5,] 9.763668 12.67313
## [6,] 10.015440 12.72405
## [7,] 9.933469 12.59941
## [8,] 9.941194 12.51213
## [9,] 9.935688 12.50187
## [10,] 9.841084 12.31474
## [11,] 9.656931 12.27403
## [12,] 9.468776 12.23455
## [13,] 9.394844 12.22656
## [14,] 9.386473 12.34728
## [15,] 9.385739 12.36164
## [16,] 9.380952 12.45537
## [17,] 9.357978 12.46236
## [18,] 9.359019 12.56527
## [19,] 9.386005 12.59624
## attr(,"bbox")
## [1] 9.357978 12.226560 10.015440 12.724050
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 19
##
##
## [[23]]
## [,1] [,2]
## [1,] 10.66777 12.64803
## [2,] 10.77667 12.65201
## [3,] 10.88756 12.64401
## [4,] 10.88849 12.49214
## [5,] 10.88444 12.34927
## [6,] 10.88438 12.19341
## [7,] 10.87632 12.03754
## [8,] 10.64090 12.01312
## [9,] 10.42425 11.99064
## [10,] 10.42175 12.44026
## [11,] 10.66970 12.48318
## [12,] 10.66734 12.58011
## [13,] 10.66697 12.59503
## [14,] 10.66777 12.64803
## attr(,"bbox")
## [1] 10.42175 11.99064 10.88849 12.65201
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 14
##
##
## [[24]]
## [,1] [,2]
## [1,] 8.148235 12.64535
## [2,] 8.223163 12.64534
## [3,] 8.296090 12.64534
## [4,] 8.299444 12.60978
## [5,] 8.301018 12.59309
## [6,] 8.305043 12.55042
## [7,] 8.332018 12.54742
## [8,] 8.348999 12.54742
## [9,] 8.347977 12.49346
## [10,] 8.350509 12.48193
## [11,] 8.358474 12.44566
## [12,] 8.369997 12.42818
## [13,] 8.391569 12.40738
## [14,] 8.406447 12.39673
## [15,] 8.438839 12.37356
## [16,] 8.481381 12.36291
## [17,] 8.502639 12.35708
## [18,] 8.523073 12.35531
## [19,] 8.554715 12.35533
## [20,] 8.565508 12.35533
## [21,] 8.589687 12.36155
## [22,] 8.605762 12.36794
## [23,] 8.657633 12.39151
## [24,] 8.694607 12.41749
## [25,] 8.709477 12.42731
## [26,] 8.708920 12.39049
## [27,] 8.708316 12.35050
## [28,] 8.707547 12.30259
## [29,] 8.767486 12.29558
## [30,] 8.758468 12.23064
## [31,] 8.784432 12.20666
## [32,] 8.790393 12.12573
## [33,] 8.643800 12.10209
## [34,] 8.538624 12.08879
## [35,] 8.519205 12.11523
## [36,] 8.505514 12.12821
## [37,] 8.495033 12.13777
## [38,] 8.471713 12.14575
## [39,] 8.434194 12.13921
## [40,] 8.408770 12.13477
## [41,] 8.358864 12.14595
## [42,] 8.326121 12.15328
## [43,] 8.310225 12.15693
## [44,] 8.250281 12.19443
## [45,] 8.225327 12.20686
## [46,] 8.190520 12.23764
## [47,] 8.174397 12.25043
## [48,] 8.157790 12.28122
## [49,] 8.148350 12.29872
## [50,] 8.141801 12.32647
## [51,] 8.136807 12.34855
## [52,] 8.131140 12.37360
## [53,] 8.137152 12.41656
## [54,] 8.149159 12.46251
## [55,] 8.160175 12.52846
## [56,] 8.155201 12.57941
## [57,] 8.147226 12.61938
## [58,] 8.148235 12.64535
## attr(,"bbox")
## [1] 8.131140 12.088790 8.790393 12.645350
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 58
##
##
## [[25]]
## [,1] [,2]
## [1,] 8.785566 12.53237
## [2,] 8.855528 12.60630
## [3,] 8.913297 12.60980
## [4,] 9.124277 12.63424
## [5,] 9.124982 12.58726
## [6,] 9.127170 12.44146
## [7,] 9.129522 12.28463
## [8,] 9.131060 12.18216
## [9,] 9.116095 12.17965
## [10,] 8.955190 12.15301
## [11,] 8.790393 12.12573
## [12,] 8.784432 12.20666
## [13,] 8.758468 12.23064
## [14,] 8.767486 12.29558
## [15,] 8.707547 12.30259
## [16,] 8.708316 12.35050
## [17,] 8.708920 12.39049
## [18,] 8.709477 12.42731
## [19,] 8.726524 12.45521
## [20,] 8.758386 12.49685
## [21,] 8.785566 12.53237
## attr(,"bbox")
## [1] 8.707547 12.125730 9.131060 12.634240
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 21
##
##
## [[26]]
## [,1] [,2]
## [1,] 9.124277 12.63424
## [2,] 9.258048 12.63061
## [3,] 9.383021 12.62721
## [4,] 9.386005 12.59624
## [5,] 9.359019 12.56527
## [6,] 9.357978 12.46236
## [7,] 9.380952 12.45537
## [8,] 9.385739 12.36164
## [9,] 9.386473 12.34728
## [10,] 9.394844 12.22656
## [11,] 9.258298 12.20358
## [12,] 9.131060 12.18216
## [13,] 9.129522 12.28463
## [14,] 9.127170 12.44146
## [15,] 9.124982 12.58726
## [16,] 9.124277 12.63424
## attr(,"bbox")
## [1] 9.124277 12.182160 9.394844 12.634240
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 16
##
##
## [[27]]
## [,1] [,2]
## [1,] 9.468776 12.23455
## [2,] 9.656931 12.27403
## [3,] 9.841084 12.31474
## [4,] 9.966332 12.34139
## [5,] 10.001230 12.17054
## [6,] 9.884351 12.19153
## [7,] 9.892879 12.17385
## [8,] 9.899952 12.16369
## [9,] 9.914780 12.14239
## [10,] 9.956293 12.08275
## [11,] 9.991961 12.03151
## [12,] 10.007160 12.00968
## [13,] 9.755397 12.00493
## [14,] 9.480667 12.00275
## [15,] 9.481962 12.06578
## [16,] 9.478416 12.16182
## [17,] 9.468776 12.23455
## attr(,"bbox")
## [1] 9.468776 12.002750 10.007160 12.341390
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 17
##
##
## [[28]]
## [,1] [,2]
## [1,] 9.394844 12.22656
## [2,] 9.468776 12.23455
## [3,] 9.478416 12.16182
## [4,] 9.481962 12.06578
## [5,] 9.480667 12.00275
## [6,] 9.482347 11.97569
## [7,] 9.488198 11.88145
## [8,] 9.489535 11.85992
## [9,] 9.492552 11.75197
## [10,] 9.343692 11.73500
## [11,] 9.200423 11.76711
## [12,] 9.084968 11.79299
## [13,] 9.090070 11.87397
## [14,] 9.093653 11.93083
## [15,] 9.098085 12.02578
## [16,] 9.100104 12.06358
## [17,] 9.103628 12.12844
## [18,] 9.115089 12.15866
## [19,] 9.131060 12.18216
## [20,] 9.258298 12.20358
## [21,] 9.394844 12.22656
## attr(,"bbox")
## [1] 9.084968 11.735000 9.492552 12.234550
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 21
##
##
## [[29]]
## [,1] [,2]
## [1,] 9.116095 12.17965
## [2,] 9.131060 12.18216
## [3,] 9.115089 12.15866
## [4,] 9.103628 12.12844
## [5,] 9.100104 12.06358
## [6,] 9.098085 12.02578
## [7,] 9.093653 11.93083
## [8,] 9.090070 11.87397
## [9,] 9.084968 11.79299
## [10,] 8.898137 11.80248
## [11,] 8.868270 12.01382
## [12,] 8.811334 12.03081
## [13,] 8.813672 12.06682
## [14,] 8.815355 12.09276
## [15,] 8.790384 12.10575
## [16,] 8.790393 12.12573
## [17,] 8.955190 12.15301
## [18,] 9.116095 12.17965
## attr(,"bbox")
## [1] 8.790384 11.792990 9.131060 12.182160
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 18
##
##
## [[30]]
## [,1] [,2]
## [1,] 8.643800 12.10209
## [2,] 8.790393 12.12573
## [3,] 8.790384 12.10575
## [4,] 8.815355 12.09276
## [5,] 8.813672 12.06682
## [6,] 8.811334 12.03081
## [7,] 8.868270 12.01382
## [8,] 8.898137 11.80248
## [9,] 8.829165 11.80486
## [10,] 8.713338 11.81002
## [11,] 8.618428 11.80504
## [12,] 8.508537 11.80904
## [13,] 8.386433 11.78549
## [14,] 8.404299 11.80047
## [15,] 8.412508 11.80735
## [16,] 8.455998 11.84761
## [17,] 8.471678 11.86373
## [18,] 8.517660 11.91471
## [19,] 8.525978 11.92466
## [20,] 8.563573 12.02385
## [21,] 8.544462 12.07360
## [22,] 8.538624 12.08879
## [23,] 8.643800 12.10209
## attr(,"bbox")
## [1] 8.386433 11.785490 8.898137 12.125730
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 23
##
##
## [[31]]
## [,1] [,2]
## [1,] 6.763314 11.98112
## [2,] 6.942171 12.05603
## [3,] 7.004114 12.06701
## [4,] 7.078046 12.07699
## [5,] 7.125001 12.07898
## [6,] 7.185831 11.84824
## [7,] 7.013006 11.82622
## [8,] 6.741261 11.79928
## [9,] 6.456532 11.78133
## [10,] 6.674380 11.93117
## [11,] 6.763314 11.98112
## attr(,"bbox")
## [1] 6.456532 11.781330 7.185831 12.078980
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 11
##
##
## [[32]]
## [,1] [,2]
## [1,] 10.42425 11.99064
## [2,] 10.64090 12.01312
## [3,] 10.87632 12.03754
## [4,] 10.88627 11.94662
## [5,] 10.91524 11.94162
## [6,] 10.93222 11.93862
## [7,] 10.96319 11.94062
## [8,] 10.98317 11.93662
## [9,] 10.99516 11.93762
## [10,] 11.04508 11.85368
## [11,] 11.11297 11.76176
## [12,] 11.20483 11.63885
## [13,] 11.12691 11.63387
## [14,] 11.07397 11.65785
## [15,] 10.97407 11.66186
## [16,] 10.95212 11.72781
## [17,] 10.87822 11.72326
## [18,] 10.80626 11.71883
## [19,] 10.78336 11.71955
## [20,] 10.76621 11.72008
## [21,] 10.67738 11.72285
## [22,] 10.65640 11.72285
## [23,] 10.61346 11.72330
## [24,] 10.56050 11.72386
## [25,] 10.53153 11.74385
## [26,] 10.42564 11.74186
## [27,] 10.42425 11.99064
## attr(,"bbox")
## [1] 10.42425 11.63387 11.20483 12.03754
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 27
##
##
## [[33]]
## [,1] [,2]
## [1,] 9.480667 12.00275
## [2,] 9.755397 12.00493
## [3,] 10.007160 12.00968
## [4,] 10.024220 11.96231
## [5,] 10.042080 11.91276
## [6,] 10.048910 11.85880
## [7,] 10.049290 11.82283
## [8,] 10.049750 11.77896
## [9,] 10.049910 11.76295
## [10,] 10.026480 11.78770
## [11,] 9.986115 11.80286
## [12,] 9.973393 11.80906
## [13,] 9.925159 11.82385
## [14,] 9.877091 11.84418
## [15,] 9.859212 11.85174
## [16,] 9.831771 11.85882
## [17,] 9.800300 11.86683
## [18,] 9.639468 11.86325
## [19,] 9.489535 11.85992
## [20,] 9.488198 11.88145
## [21,] 9.482347 11.97569
## [22,] 9.480667 12.00275
## attr(,"bbox")
## [1] 9.480667 11.762950 10.049910 12.009680
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 22
##
##
## [[34]]
## [,1] [,2]
## [1,] 7.185831 11.84824
## [2,] 7.317721 11.85615
## [3,] 7.508541 11.87212
## [4,] 7.484550 11.83815
## [5,] 7.463559 11.80918
## [6,] 7.469532 11.75722
## [7,] 7.493496 11.72825
## [8,] 7.492452 11.61934
## [9,] 7.688279 11.66328
## [10,] 7.733179 11.52739
## [11,] 7.710183 11.48244
## [12,] 7.674203 11.44647
## [13,] 7.624242 11.42150
## [14,] 7.562289 11.38854
## [15,] 7.468364 11.34958
## [16,] 7.398374 11.32995
## [17,] 7.375505 11.47149
## [18,] 7.315612 11.58639
## [19,] 7.123811 11.53079
## [20,] 6.984895 11.49052
## [21,] 6.971885 11.63357
## [22,] 6.966994 11.68735
## [23,] 7.019946 11.69534
## [24,] 7.013006 11.82622
## [25,] 7.185831 11.84824
## attr(,"bbox")
## [1] 6.966994 11.329950 7.733179 11.872120
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 25
##
##
## [[35]]
## [,1] [,2]
## [1,] 9.489535 11.85992
## [2,] 9.639468 11.86325
## [3,] 9.800300 11.86683
## [4,] 9.831771 11.85882
## [5,] 9.859212 11.85174
## [6,] 9.877091 11.84418
## [7,] 9.925159 11.82385
## [8,] 9.973393 11.80906
## [9,] 9.986115 11.80286
## [10,] 10.026480 11.78770
## [11,] 10.049910 11.76295
## [12,] 10.056590 11.74391
## [13,] 9.917131 11.73793
## [14,] 9.918046 11.53111
## [15,] 9.500506 11.54202
## [16,] 9.341615 11.54617
## [17,] 9.343692 11.73500
## [18,] 9.492552 11.75197
## [19,] 9.489535 11.85992
## attr(,"bbox")
## [1] 9.341615 11.531110 10.056590 11.866830
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 19
##
##
## [[36]]
## [,1] [,2]
## [1,] 6.741261 11.79928
## [2,] 7.013006 11.82622
## [3,] 7.019946 11.69534
## [4,] 6.966994 11.68735
## [5,] 6.971885 11.63357
## [6,] 6.984895 11.49052
## [7,] 6.891401 11.47941
## [8,] 6.777808 11.46591
## [9,] 6.758334 11.46233
## [10,] 6.678177 11.44760
## [11,] 6.610277 11.44478
## [12,] 6.558312 11.44220
## [13,] 6.481366 11.43963
## [14,] 6.456532 11.78133
## [15,] 6.741261 11.79928
## attr(,"bbox")
## [1] 6.456532 11.439630 7.019946 11.826220
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 15
##
##
## [[37]]
## [,1] [,2]
## [1,] 8.713338 11.81002
## [2,] 8.829165 11.80486
## [3,] 8.898137 11.80248
## [4,] 9.084968 11.79299
## [5,] 9.085843 11.63908
## [6,] 9.095860 11.55919
## [7,] 9.082836 11.45633
## [8,] 9.000834 11.45802
## [9,] 8.935976 11.46030
## [10,] 8.922977 11.43432
## [11,] 8.772132 11.45033
## [12,] 8.713338 11.81002
## attr(,"bbox")
## [1] 8.713338 11.434320 9.095860 11.810020
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 12
##
##
## [[38]]
## [,1] [,2]
## [1,] 9.084968 11.79299
## [2,] 9.200423 11.76711
## [3,] 9.343692 11.73500
## [4,] 9.341615 11.54617
## [5,] 9.207753 11.55326
## [6,] 9.095860 11.55919
## [7,] 9.085843 11.63908
## [8,] 9.084968 11.79299
## attr(,"bbox")
## [1] 9.084968 11.546170 9.343692 11.792990
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 8
##
##
## [[39]]
## [,1] [,2]
## [1,] 6.181706 11.55357
## [2,] 6.394562 11.71040
## [3,] 6.456532 11.78133
## [4,] 6.481366 11.43963
## [5,] 6.336500 11.42066
## [6,] 6.335865 11.39568
## [7,] 6.335484 11.38070
## [8,] 6.356451 11.34972
## [9,] 6.359416 11.27379
## [10,] 6.316442 11.23283
## [11,] 6.167586 11.23085
## [12,] 6.179293 11.18159
## [13,] 6.180611 11.16772
## [14,] 6.212476 11.12212
## [15,] 6.220480 11.09996
## [16,] 6.220463 11.05700
## [17,] 6.088609 11.09798
## [18,] 6.088623 11.13295
## [19,] 6.084633 11.14894
## [20,] 6.085638 11.16092
## [21,] 6.035704 11.20289
## [22,] 5.984754 11.20390
## [23,] 5.948804 11.23987
## [24,] 5.969793 11.26385
## [25,] 5.970815 11.31780
## [26,] 5.874907 11.31382
## [27,] 6.181706 11.55357
## attr(,"bbox")
## [1] 5.874907 11.057000 6.481366 11.781330
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 27
##
##
## [[40]]
## [,1] [,2]
## [1,] 10.05659 11.74391
## [2,] 10.25280 11.74282
## [3,] 10.42564 11.74186
## [4,] 10.53153 11.74385
## [5,] 10.56050 11.72386
## [6,] 10.61346 11.72330
## [7,] 10.65640 11.72285
## [8,] 10.67738 11.72285
## [9,] 10.76621 11.72008
## [10,] 10.78336 11.71955
## [11,] 10.80626 11.71883
## [12,] 10.80224 11.67487
## [13,] 10.80322 11.66159
## [14,] 10.80609 11.62091
## [15,] 10.80820 11.59095
## [16,] 10.80869 11.57467
## [17,] 10.80493 11.55899
## [18,] 10.79919 11.53500
## [19,] 10.80214 11.43408
## [20,] 10.80206 11.23226
## [21,] 10.73412 11.22827
## [22,] 10.68717 11.22229
## [23,] 10.61324 11.21630
## [24,] 10.59525 11.21331
## [25,] 10.55529 11.22030
## [26,] 10.48037 11.22131
## [27,] 10.42143 11.23231
## [28,] 10.36449 11.25430
## [29,] 10.34352 11.27228
## [30,] 10.32956 11.32824
## [31,] 10.32059 11.36920
## [32,] 10.27964 11.40318
## [33,] 10.22671 11.45614
## [34,] 10.16880 11.52309
## [35,] 10.10988 11.59003
## [36,] 10.07095 11.67197
## [37,] 10.05659 11.74391
## attr(,"bbox")
## [1] 10.05659 11.21331 10.80869 11.74391
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 37
##
##
## [[41]]
## [,1] [,2]
## [1,] 10.80626 11.71883
## [2,] 10.87822 11.72326
## [3,] 10.95212 11.72781
## [4,] 10.97407 11.66186
## [5,] 11.07397 11.65785
## [6,] 11.12691 11.63387
## [7,] 11.12689 11.59590
## [8,] 11.13887 11.56992
## [9,] 11.13985 11.53096
## [10,] 11.14184 11.50598
## [11,] 11.12285 11.49299
## [12,] 11.12279 11.33613
## [13,] 11.12221 11.31588
## [14,] 11.09847 11.30317
## [15,] 11.07896 11.29528
## [16,] 11.06859 11.29187
## [17,] 11.01487 11.27420
## [18,] 10.96391 11.26521
## [19,] 10.91995 11.25423
## [20,] 10.90186 11.25213
## [21,] 10.87750 11.24687
## [22,] 10.84074 11.23894
## [23,] 10.81888 11.23517
## [24,] 10.80206 11.23226
## [25,] 10.80214 11.43408
## [26,] 10.79919 11.53500
## [27,] 10.80493 11.55899
## [28,] 10.80869 11.57467
## [29,] 10.80820 11.59095
## [30,] 10.80609 11.62091
## [31,] 10.80322 11.66159
## [32,] 10.80224 11.67487
## [33,] 10.80626 11.71883
## attr(,"bbox")
## [1] 10.79919 11.23226 11.14184 11.72781
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 33
##
##
## [[42]]
## [,1] [,2]
## [1,] 6.758334 11.46233
## [2,] 6.777808 11.46591
## [3,] 6.891401 11.47941
## [4,] 6.984895 11.49052
## [5,] 7.123811 11.53079
## [6,] 7.315612 11.58639
## [7,] 7.375505 11.47149
## [8,] 7.398374 11.32995
## [9,] 7.404403 11.29264
## [10,] 7.382380 11.18474
## [11,] 7.359378 11.12479
## [12,] 7.327387 11.07185
## [13,] 7.275407 10.99792
## [14,] 6.827894 11.11687
## [15,] 6.879882 11.21178
## [16,] 6.811956 11.22877
## [17,] 6.830951 11.26274
## [18,] 6.771011 11.26575
## [19,] 6.754039 11.29272
## [20,] 6.786021 11.32569
## [21,] 6.778052 11.38064
## [22,] 6.758076 11.39263
## [23,] 6.758334 11.46233
## attr(,"bbox")
## [1] 6.754039 10.997920 7.404403 11.586390
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 23
##
##
## [[43]]
## [,1] [,2]
## [1,] 9.095860 11.55919
## [2,] 9.207753 11.55326
## [3,] 9.341615 11.54617
## [4,] 9.338522 11.31238
## [5,] 9.207670 11.31017
## [6,] 9.103749 11.30841
## [7,] 9.082836 11.45633
## [8,] 9.095860 11.55919
## attr(,"bbox")
## [1] 9.082836 11.308410 9.341615 11.559190
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 8
##
##
## [[44]]
## [,1] [,2]
## [1,] 9.341615 11.54617
## [2,] 9.500506 11.54202
## [3,] 9.918046 11.53111
## [4,] 9.963892 11.26834
## [5,] 9.778072 11.26536
## [6,] 9.776052 11.21141
## [7,] 9.656154 11.21248
## [8,] 9.548276 11.21344
## [9,] 9.546307 11.28537
## [10,] 9.391455 11.27540
## [11,] 9.335508 11.27241
## [12,] 9.338522 11.31238
## [13,] 9.341615 11.54617
## attr(,"bbox")
## [1] 9.335508 11.211410 9.963892 11.546170
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 13
##
##
## [[45]]
## [,1] [,2]
## [1,] 9.000834 11.45802
## [2,] 9.082836 11.45633
## [3,] 9.103749 11.30841
## [4,] 9.111679 11.15754
## [5,] 9.107675 11.13556
## [6,] 9.098839 11.11539
## [7,] 9.093675 11.10359
## [8,] 9.171593 11.09159
## [9,] 9.174187 11.06031
## [10,] 9.175572 11.04863
## [11,] 9.178549 10.99967
## [12,] 9.146570 10.97670
## [13,] 9.135557 10.91775
## [14,] 9.100430 10.89909
## [15,] 9.095559 10.82883
## [16,] 8.850796 10.82187
## [17,] 8.815848 10.86784
## [18,] 8.789886 10.89981
## [19,] 8.765927 10.93878
## [20,] 8.752955 10.97875
## [21,] 8.747973 11.01072
## [22,] 8.753979 11.03570
## [23,] 8.750992 11.06067
## [24,] 8.622117 11.05969
## [25,] 8.631121 11.08966
## [26,] 8.634133 11.12663
## [27,] 8.632165 11.19957
## [28,] 8.665152 11.24652
## [29,] 8.683140 11.25751
## [30,] 8.697130 11.26750
## [31,] 8.675208 11.40538
## [32,] 8.772132 11.45033
## [33,] 8.922977 11.43432
## [34,] 8.935976 11.46030
## [35,] 9.000834 11.45802
## attr(,"bbox")
## [1] 8.622117 10.821870 9.178549 11.460300
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 35
##
##
## [[46]]
## [,1] [,2]
## [1,] 6.481366 11.43963
## [2,] 6.558312 11.44220
## [3,] 6.610277 11.44478
## [4,] 6.678177 11.44760
## [5,] 6.649185 11.39964
## [6,] 6.647165 11.34469
## [7,] 6.639161 11.31572
## [8,] 6.522276 11.31873
## [9,] 6.520258 11.26977
## [10,] 6.516248 11.23680
## [11,] 6.542212 11.20982
## [12,] 6.549193 11.18085
## [13,] 6.555157 11.10791
## [14,] 6.534173 11.09692
## [15,] 6.537121 10.97803
## [16,] 6.220463 11.05700
## [17,] 6.220480 11.09996
## [18,] 6.212476 11.12212
## [19,] 6.180611 11.16772
## [20,] 6.179293 11.18159
## [21,] 6.167586 11.23085
## [22,] 6.316442 11.23283
## [23,] 6.359416 11.27379
## [24,] 6.356451 11.34972
## [25,] 6.335484 11.38070
## [26,] 6.335865 11.39568
## [27,] 6.336500 11.42066
## [28,] 6.481366 11.43963
## attr(,"bbox")
## [1] 6.167586 10.978030 6.678177 11.447600
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 28
##
##
## [[47]]
## [,1] [,2]
## [1,] 11.09847 11.30317
## [2,] 11.12221 11.31588
## [3,] 11.12076 11.26519
## [4,] 11.11973 11.19925
## [5,] 11.12370 11.12731
## [6,] 11.14465 11.05637
## [7,] 11.18160 11.02540
## [8,] 11.22154 10.97843
## [9,] 11.25649 10.92248
## [10,] 11.27545 10.88850
## [11,] 11.28742 10.84254
## [12,] 11.28640 10.79059
## [13,] 11.01367 10.78863
## [14,] 10.96372 10.79962
## [15,] 10.88380 10.80762
## [16,] 10.78190 10.81263
## [17,] 10.70898 10.81964
## [18,] 10.69400 10.83363
## [19,] 10.64907 10.89058
## [20,] 10.61611 10.92556
## [21,] 10.58816 10.97851
## [22,] 10.59217 10.99450
## [23,] 10.69113 11.14036
## [24,] 10.68717 11.22229
## [25,] 10.73412 11.22827
## [26,] 10.80206 11.23226
## [27,] 10.81888 11.23517
## [28,] 10.84074 11.23894
## [29,] 10.87750 11.24687
## [30,] 10.90186 11.25213
## [31,] 10.91995 11.25423
## [32,] 10.96391 11.26521
## [33,] 11.01487 11.27420
## [34,] 11.06859 11.29187
## [35,] 11.07896 11.29528
## [36,] 11.09847 11.30317
## attr(,"bbox")
## [1] 10.58816 10.78863 11.28742 11.31588
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 36
##
##
## [[48]]
## [,1] [,2]
## [1,] 9.207670 11.31017
## [2,] 9.338522 11.31238
## [3,] 9.335508 11.27241
## [4,] 9.391455 11.27540
## [5,] 9.393032 11.26274
## [6,] 9.393250 11.25140
## [7,] 9.394897 11.16547
## [8,] 9.394545 11.13953
## [9,] 9.393360 11.05260
## [10,] 9.282284 11.04803
## [11,] 9.175572 11.04863
## [12,] 9.174187 11.06031
## [13,] 9.171593 11.09159
## [14,] 9.093675 11.10359
## [15,] 9.098839 11.11539
## [16,] 9.107675 11.13556
## [17,] 9.111679 11.15754
## [18,] 9.103749 11.30841
## [19,] 9.207670 11.31017
## attr(,"bbox")
## [1] 9.093675 11.048030 9.394897 11.312380
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 19
##
##
## [[49]]
## [,1] [,2]
## [1,] 9.391455 11.27540
## [2,] 9.546307 11.28537
## [3,] 9.548276 11.21344
## [4,] 9.656154 11.21248
## [5,] 9.776052 11.21141
## [6,] 9.779012 11.11849
## [7,] 9.774969 11.00759
## [8,] 9.775935 10.92866
## [9,] 9.781898 10.85373
## [10,] 9.095559 10.82883
## [11,] 9.100430 10.89909
## [12,] 9.135557 10.91775
## [13,] 9.146570 10.97670
## [14,] 9.178549 10.99967
## [15,] 9.175572 11.04863
## [16,] 9.282284 11.04803
## [17,] 9.393360 11.05260
## [18,] 9.394545 11.13953
## [19,] 9.394897 11.16547
## [20,] 9.393250 11.25140
## [21,] 9.393032 11.26274
## [22,] 9.391455 11.27540
## attr(,"bbox")
## [1] 9.095559 10.828830 9.781898 11.285370
## attr(,"ringDir")
## [1] 1
## attr(,"after")
## [1] NA
## attr(,"plotOrder")
## [1] 1
## attr(,"nParts")
## [1] 1
## attr(,"pstart")
## attr(,"pstart")$from
## [1] 1
##
## attr(,"pstart")$to
## [1] 22
##
##
## attr(,"region.id")
## [1] "1005" "1001" "1006" "1002" "1007" "1008" "1004" "1003" "1018" "1010"
## [11] "1038" "1037" "1039" "1040" "1009" "1036" "1011" "1042" "1041" "1017"
## [21] "1043" "1019" "1012" "1035" "1032" "1020" "1021" "1031" "1033" "1034"
## [31] "1045" "1013" "1022" "1044" "1023" "1046" "1030" "1024" "1047" "1016"
## [41] "1014" "1049" "1029" "1025" "1028" "1048" "1015" "1027" "1026"
## attr(,"class")
## [1] "polylist"
## attr(,"maplim")
## attr(,"maplim")$x
## [1] 5.874907 11.287420
##
## attr(,"maplim")$y
## [1] 10.78863 14.74245
##
## attr(,"plotOrder")
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## [26] 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## attr(,"after")
## [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## [26] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NACada elemento de la lista contiene una matriz de dos columnas con coordenadas x en la columna 1 y coordenadas y en la columna 2, con tantas filas como puntos digitalizados haya en el contorno del polígono en cuestión. Después de la matriz de coordenadas viene el cuadro de límite y varias opciones de trazado.
Hay una función muy útil \(poly2nb\) que toma la lista de polígonos y determina qué regiones son vecinas entre sí buscando límites compartidos. El resultado es un objeto nb (aquí llamado colnbs), y podemos obtener una verificación visual de qué tan bien ha funcionado poly2nb superponiendo las relaciones vecinas en un mapa de los contornos del polígono:
#colnbs <- poly2nb(polys)
#plot(c(5.5,11.5),c(10.5,15),type="n",xlab="",ylab="")
#for (i in 1:49) polygon(polys[[i]][,1],polys[[i]][,2],col="lightgrey")
#plot(colnbs,coords,add=T,col="red")26.5 - Información geoestadística
Los datos mapeados comúnmente muestran el valor de una variable de respuesta continua (por ejemplo, la concentración de un mineral) en diferentes ubicaciones espaciales. El problema fundamental con este tipo de datos es la pseudorreplicación espacial. Los puntos calientes tienden a generar una gran cantidad de datos, y estos datos tienden a ser bastante similares porque provienen esencialmente del mismo lugar. Los puntos fríos están poco representados y, por lo general, están muy separados. Grandes áreas entre los puntos fríos no tienen ningún dato.
La estadística espacial tiene en cuenta esta autocorrelación espacial de diversas formas. La herramienta fundamental de la estadística espacial es el variograma (o semivariograma). Esto mide la rapidez con la que la autocorrelación espacial, γ (h), cae al aumentar la distancia:
\[\gamma(h)=\frac{1}{2|N(h)|}\sum_{N(h)}{(z_i-z_j)^2}\]
Aquí N (h) es el conjunto de todas las distancias euclidianas por pares i - j = h, | N (h) | es el número de pares distintos dentro de N (h), y zi y zj son valores de la variable de respuesta en las ubicaciones espaciales i y j. Hay dos reglas generales importantes: (1) la distancia de confiabilidad del variograma es menos de la mitad de la distancia máxima en todo el campo de datos; y (2) solo debe considerar producir un variograma empírico cuando tenga más de 30 puntos de datos en el mapa.
Las gráficas del variograma empírico contra la distancia se caracterizan por algunas características con nombres curiosos que delatan su origen en la prospección geológica:
nugget variación a pequeña escala más error de medición
sill el valor asintótico de \(\gamma(h)\) como \(h\rightarrow \infty\), que representa la varianza del campo aleatorio
range la distancia umbral (si existe) más allá de la cual los datos ya no están autocorrelacionados.
Las gráficas de variograma que no tienen asíntotas pueden ser sintomáticas de datos con tendencia o un proceso estocástico no estacionario. El covariograma C (h) es la covarianza de los valores de z en la separación h, para todo \(i\) e \(i + h\) dentro de la distancia máxima en todo el campo de datos:
\[cov(Z(i+h),Z(i))=C(h)\]
El correlograma es una relación de covarianzas:
\[\rho(h)=\frac{C(h)}{C(0)}=1-\frac{\gamma(h)}{C(0)}\]
Aquí \(C(0)\) es la varianza del campo aleatorio y \(\gamma(h)\) es el variograma. Donde el variograma aumenta con la distancia, el correlograma y el covariograma disminuyen con la distancia.
El variograma asume que los datos están sin tendencia. Si hay tendencias, entonces una opción es el pulido medio. Esto implica modelar efectos de fila y columna del mapa de esta manera:
\[y \sim overall~mean+row~effect+column~effect+residual\]
Este modelo bidireccional asume efectos aditivos y no funcionaría si hubiera una interacción entre las filas y columnas del mapa. Una alternativa sería utilizar un modelo aditivo generalizado con suavizadores no paramétricos para latitud y longitud.
La anisotropía ocurre cuando la autocorrelación espacial cambia con la dirección. Si el umbral cambia con la dirección, esto se llama anisotropía zonal. Cuando es el rango el que cambia con la dirección, el proceso se llama anisotropía geométrica.
Los geógrafos tienen la maravillosa habilidad de hacer que las ideas más simples parezcan complicadas. Kriging no es más que una interpolación lineal a través del espacio. El kriging ordinario utiliza un modelo de función aleatoria de correlación espacial para calcular una combinación lineal ponderada de las muestras disponibles para predecir la respuesta para una ubicación no medida. El kriging universal es una modificación del kriging ordinario que permite las tendencias espaciales. No hablamos más aquí de los modelos de predicción espacial; los detalles se pueden encontrar en Kaluzny et al. (1998). Nuestra preocupación es el uso de información espacial en la interpretación de estudios experimentales u observacionales que tienen una única variable de respuesta. El énfasis está en el uso de medidas específicas de la ubicación para modelar la estructura de autocorrelación espacial de los datos.
La idea de un variograma es ilustrar la forma en que la varianza espacial aumenta con la escala espacial (o, alternativamente, cómo la correlación entre vecinos disminuye con la distancia). Confusamente, R tiene dos funciones con el mismo nombre: variograma (“v” minúscula) está en la biblioteca espacial y Variograma (“V” mayúscula) está en nlme. Su uso se contrasta aquí para los datos de \(ragwort\).
Para usar \(variogram\) de la biblioteca \(Spatial\), necesita crear una superficie de tendencia o un objeto kriging con columnas $ x $, $ y $ y $ z $. Las dos primeras columnas son las coordenadas espaciales, mientras que la tercera contiene la variable de respuesta (diámetro del tallo basal en el caso de los datos de \(ragwort\))
names(data_rwm)## [1] "ROW" "COLUMN" "X" "Y" "rosette"
## [6] "regenerating" "stems" "diameter" "xcoord" "ycoord"
## [11] "type"diameter_rgw <- as.factor(data_rwm$diameter)
dts <- data.frame(x=xcoord,y=ycoord,z=diameter_rgw)A continuación, debe crear una superficie de tendencia utilizando una función como \(surf.ls\):
#surface <- surf.ls(2,dts)Este objeto de superficie de tendencia es entonces el primer argumento del variograma, seguido del número de bins (aquí 300). La función calcula la diferencia cuadrática promedio para pares con separación en cada contenedor, y devuelve resultados para contenedores que contienen seis o más pares:
#variogram(surface,300)La función hermana es el correlograma, que toma argumentos idénticos:
#correlogram(surface,300)Las correlaciones positivas han desaparecido en unos 100 cm. Las correlaciones en xp = 3000 son efectos de borde espurios.
Para la función Variograma en la biblioteca nlme, necesita ajustar un modelo (normalmente usando gls o lme), luego proporcionar el objeto modelo junto con una función de formulario en la llamada:
#model_37 <- gls(diameter~xcoord+ycoord)
#plot(Variogram(model_37,form= ~xcoord+ycoord))26.6 - Modelos de regresión con errores espacialmente correlacionados: mínimos cuadrados generalizados
Analizamos el uso de modelos lineales de efectos mixtos para tratar los efectos aleatorios y la pseudorreplicación temporal. Aquí ilustramos el uso de mínimos cuadrados generalizados (GLS) para modelos de regresión donde esperaríamos que los valores vecinos de la variable de respuesta estuvieran correlacionados. La gran ventaja de la función gls es que los errores pueden correlacionarse y / o tener variaciones desiguales. La función gls es parte del paquete \(nlme\):
El siguiente ejemplo es un ensayo a escala geográfica para comparar los rendimientos de 56 variedades diferentes de trigo. Lo que hace que el análisis sea más desafiante es que las granjas que llevaron a cabo el ensayo se distribuyeron en una amplia gama de latitudes y longitudes.
spatialdata <- read.delim("D:/Kevin/Trabajos/Diseno de Experimentos/therbook/spatialdata.txt",header=T)
attach(spatialdata)
names(spatialdata)## [1] "Block" "variety" "yield" "latitude" "longitude"Comenzamos con la inspección de datos gráficos para ver el efecto de la ubicación en el rendimiento:
fig_sp <- plot_ly(data=data_sp,x = latitude_sp,y = yield_sp)
fig_sp <- fig_sp %>% add_markers(color ="black")
fig_sp <- fig_sp %>% layout(xaxis=list(title="Latitude"),
yaxis = list(title="Yield"))
fig_sp## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsfig_sp_l <- plot_ly(data=data_sp,x = longitude_sp,y = yield_sp)
fig_sp_l <- fig_sp_l %>% add_markers(color ="black")
fig_sp_l <- fig_sp_l %>% layout(xaxis=list(title="Longitude"),
yaxis = list(title="Yield"))
fig_sp_l## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levelsClaramente, existen grandes efectos de la latitud y la longitud tanto en el rendimiento medio como en la variación del rendimiento. El efecto de latitud parece un efecto de umbral, con poco impacto para latitudes inferiores a 30. El efecto de longitud parece más continuo pero hay un indicio de no linealidad (tal vez incluso una joroba). Las variedades difieren sustancialmente en sus rendimientos medios:
yield_sp <- as.numeric(yield_sp)
variety_sp <- as.numeric(variety_sp)
yie_v <- sort(tapply(yield_sp,variety_sp,mean))
barplot(yie_v,col="green")
Las variedades de menor rendimiento producen alrededor de 20 y las más altas alrededor de 30 kg de grano por unidad de área. También hay efectos de bloque sustanciales sobre el rendimiento:
yield_sp <- as.numeric(yield_sp)
Block_sp <- as.numeric(Block_sp)
tapply(yield_sp,Block_sp,mean)## 1 2 3 4
## 99.39286 109.92727 86.69643 62.98246Aquí está el análisis más simple posible: un análisis de varianza unidireccional que usa la variedad como la única variable explicativa:
variety_sp <- as.numeric(variety_sp)
yield_sp <- as.numeric(yield_sp)
model_38 <- aov(yield_sp~variety_sp)
summary(model_38)## Df Sum Sq Mean Sq F value Pr(>F)
## variety_sp 1 666 665.7 0.293 0.589
## Residuals 222 505102 2275.2Esto dice que no hay diferencias significativas entre los rendimientos de las 56 variedades. Podemos probar un análisis de parcelas divididas usando variedades anidadas dentro de bloques:
Block_sp <- as.numeric(Block_sp)
class(yield_sp)## [1] "numeric"model_39 <- aov(yield_sp~Block_sp+variety_sp+Error(Block_sp))
summary(model_39)##
## Error: Block_sp
## Df Sum Sq Mean Sq
## Block_sp 1 49567 49567
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## variety_sp 1 686 685.6 0.333 0.565
## Residuals 221 455515 2061.2Esto no ha cambiado nuestra interpretación. Podríamos ajustar la latitud y la longitud como covariables:
model_40 <- aov(yield_sp~Block_sp+variety_sp+latitude_sp+longitude_sp)
summary(model_40)## Df Sum Sq Mean Sq F value Pr(>F)
## Block_sp 1 49567 49567 44.063 3.24e-10 ***
## variety_sp 1 686 686 0.610 0.436
## latitude_sp 10 78744 7874 7.000 2.55e-09 ***
## longitude_sp 21 163039 7764 6.902 1.47e-14 ***
## Residuals 190 213733 1125
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Esto marca una enorme diferencia. Ahora las diferencias entre variedades son cercanas a la significancia (\(p = 0.0565\)).
Finalmente, podríamos usar un modelo GLS para introducir la covarianza espacial entre los rendimientos de ubicaciones cercanas. Comenzamos haciendo un objeto de datos agrupado:
space <- groupedData(yield_sp~variety_sp|Block_sp)Usamos esto para ajustar un modelo usando \(gls\) que permite que los errores estén correlacionados y tengan variaciones desiguales. Agregaremos estas sofisticaciones más adelante:
model_41 <- gls(yield_sp~variety_sp-1,space)
summary(model_41)## Generalized least squares fit by REML
## Model: yield_sp ~ variety_sp - 1
## Data: space
## AIC BIC logLik
## 2518.041 2524.855 -1257.02
##
## Coefficients:
## Value Std.Error t-value p-value
## variety_sp 2.351217 0.1346576 17.4607 0
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.7516922 -0.4036080 0.3798594 1.0024310 2.5487185
##
## Residual standard error: 66.03223
## Degrees of freedom: 224 total; 223 residualy así sucesivamente, para las 56 variedades. Se dan las medias de variedad, en lugar de diferencias entre medias, porque eliminamos la intersección del modelo utilizando rendimiento ~ variedad-1 en lugar de rendimiento ~ variedad en la fórmula del modelo. Ahora queremos incluir la covarianza espacial. La función Variograma se aplica a model_41 así:
plot(Variogram(model_41,form=~latitude_sp+longitude_sp))
El variograma de muestra aumenta con la distancia, lo que ilustra la correlación espacial esperada. Extrapolando de nuevo a la distancia cero, parece haber una pepita de alrededor de 0,2. Hay varias suposiciones que podríamos hacer sobre la correlación espacial en estos datos. Por ejemplo, podríamos probar una estructura de correlación esférica, usando la clase corSpher (el rango de opciones para la estructura de correlación espacial se muestra en la Tabla 26.1). Necesitamos especificar la distancia a que el semivariograma alcanza primero 1. La inspección muestra que esta distancia es aproximadamente 28. Podemos actualizar model4 para incluir esta información:
model_42 <- update(model_41,corr=corSpher(c(28,0.2),form=~latitude_sp+longitude_sp,nugget=T))
summary(model_42)## Generalized least squares fit by REML
## Model: yield_sp ~ variety_sp - 1
## Data: space
## AIC BIC logLik
## 2297.296 2310.924 -1144.648
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~latitude_sp + longitude_sp
## Parameter estimate(s):
## range nugget
## 34.90009067 0.01312453
##
## Coefficients:
## Value Std.Error t-value p-value
## variety_sp -0.03117947 0.1504567 -0.2072322 0.836
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## 0.01074539 0.34714616 0.57002348 0.83476847 1.11395679
##
## Residual standard error: 156.8997
## Degrees of freedom: 224 total; 223 residualEsta es una gran mejora y el AIC ha caído de 1354,742 a 1185,863. El \(range\) (27,46) y la$ nugget$ (0,209) están muy cerca de nuestras estimaciones visuales. Existen otros tipos de modelos espaciales, por supuesto. Podríamos probar un modelo cuadrático racional (\(corRatio\)); esto necesita una estimación de la distancia a la que se encuentra el semivariograma \((1 + nugget) / 2 = 1,2 / 2 = 0,6\), así como una estimación de la pepita. La inspección da una distancia de aproximadamente 12,5, por lo que escribimos:
model_43 <- update(model_41,
corr=corRatio(c(12.5,0.2),form=~latitude_sp+longitude_sp,nugget=T))Podemos usar anova para comparar los dos modelos espaciales:
anova(model_42,model_43)## Model df AIC BIC logLik
## model_42 1 4 2297.296 2310.924 -1144.648
## model_43 2 4 2297.483 2311.111 -1144.741El modelo cuadrático racional (modelo_43) tiene el AIC más bajo y, por lo tanto, se prefiere al modelo esférico. Para probar la significancia de los parámetros de correlación espacial, necesitamos comparar el modelo espacial preferido6 con el modelo no espacial model_41 (que asumió errores espacialmente independientes):
anova(model_41,model_43)## Model df AIC BIC logLik Test L.Ratio p-value
## model_41 1 2 2518.041 2524.855 -1257.020
## model_43 2 4 2297.483 2311.111 -1144.741 1 vs 2 224.5581 <.0001Los dos grados de libertad extra utilizados para dar cuenta de la estructura espacial están claramente justificados. Necesitamos comprobar la idoneidad del modelo de \(corRatio\). Esto se hace mediante la inspección del variograma de muestra para los residuos normalizados del model_43:
plot(Variogram(model_43,resType="n"))## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## pseudoinverse used at 0.995## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## neighborhood radius 0.41921## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## There are other near singularities as well. 0.074501## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## pseudoinverse used at 0.995## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## neighborhood radius 0.41921## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## There are other near singularities as well. 0.074501## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## pseudoinverse used at 0.995## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## neighborhood radius 0.41921## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## There are other near singularities as well. 0.074501## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## pseudoinverse used at 0.995## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## neighborhood radius 0.41921## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = FALSE, :
## There are other near singularities as well. 0.074501
No hay patrón en la gráfica del variograma de muestra, por lo que concluimos que el cuadrático racional es adecuado. Para verificar la constancia de la varianza, podemos graficar los residuales normalizados contra los valores ajustados de esta manera:
plot(model_43,resid( ., type="n")~fitted(.),abline=0)
y la gráfica normal se obtiene de la forma habitual:
qqnorm(model_43,~resid(.,type="n"))
El modelo luce bien. El siguiente paso es investigar la importancia de las diferencias entre las variedades. Utilice la actualización para cambiar la estructura del modelo de \(yield~variety-1\) a \(yield~variety\):
model_44 <- update(model_43,model=yield~variety)
anova(model_44)## Denom. DF: 168
## numDF F-value p-value
## (Intercept) 1 33.90671 <.0001
## variety 55 1.28563 0.1148Las diferencias entre las variedades ahora parecen ser altamente significativas (recuerde que solo fueron marginalmente significativas con nuestro modelo lineal3 usando análisis de covarianza para tener en cuenta los efectos de latitud y longitud). Los contrastes específicos entre variedades se pueden realizar utilizando el argumento L de anova. Suponga que queremos comparar los rendimientos medios de la primera y la tercera variedad. Para hacer esto, configuramos un vector de coeficientes de contraste c (-1,0,1) y aplicamos el contraste así:
#anova(model_43,L=c(-1,0,1))Tenga en cuenta que usamos model_43 (con todas las medias de variedad), no model7 (con una intersección y contrastes de Helmert). Las variedades especificadas, Arapahoe y Buckskin, presentan diferencias muy significativas en el rendimiento medio.
26.7 - Crear un mapa de distribución de puntos a partir de una base de datos relacional
A continuación se muestra un ejemplo de cómo extraer un subconjunto relativamente pequeño de datos de una gran base de datos relacional y utilizar la información para producir un mapa de distribución de puntos. La base de datos de Access contiene dos tablas relacionadas:
\(sites\) contiene información sobre 2628 ubicaciones.
\(records\) contienen listas de especies encontradas en cada sitio (43 001 en total).
Las dos tablas están relacionadas por una variable llamada número de sitio. La tarea consiste en extraer este y norte para cada registro de una especie nombrada, y usarlos para producir un mapa de distribución de puntos, con un punto para cada sitio en el que se registró esa especie en particular.
#channel <- odbcConnect("berks")Para completar el mapa necesitas:
Leer un archivo de coordenadas \(x\) y \(y\) para el contorno de la región que se está mapeando (el condado de Berkshire en este ejemplo).
Dibuje el contorno en una parcela sin ningún rótulo en los ejes.
Utilice el eje para etiquetar con un solo dígito (las referencias de la cuadrícula tienen 5 dígitos).
Agregue una cuadrícula usando abline con líneas grises.
Especificar una especie para cartografiar.
Leer las coordenadas de esta especie en R de la base de datos de Access
Use puntos para agregar los puntos de distribución al mapa. Aquí está el esquema del condado de Berkshire.