Indices de Agregación

Objetivo

Calcular estadísticas descriptivas (media, varianza, totales) de lecturas estructuradas. Asimismo, calcular algunos índices de agregación, creación de gráficas y estimar tamaños de muestra.

Datos

• datosMuestreoCompleto.csv, cargado como dcitri

• ssdcD.R, script para calcular tamaños de muestra con la ecuación de Taylor

Los archivos se encuentran aquí:

https://sites.google.com/site/digitcognem/boxfile_pract

Paquetes

• dplyr para procesamiento de datos (group_by())

• ggplot para generar gráficas

• reshape2 para reestructurar el conjunto de datos con la función melt 

El archivo de datos contiene las capturas de Diaphorina citri durante doce meses. Se pusieron trampas amarillas en árboles de limonaria, Murraya paniculata en la ciudad de Cuitláhuac para estudiar la fluctuación poblacional de este insecto (Hernandez-Landa et al. 2013). D. citri es el vector del Huanglongbing, una enfermedad letal de los cítricos.

La estructura de las columnas es como sigue:

longitud y latitud son las coordenadas geográficas de las trampas

trampa, es el Id de la trampa

1, 2, 3,..., 12 son los meses de muestreo (Febrero 2012, ..., Enero 2013)

Total, es el total de capturas en todos los meses.

Manejo de Datos

Para procesar la información, los conteos de los doce meses exclusivamente, es necesario eliminar las columnas innecesarias: latitud, longitud, trampa yTotal. 

dcitr$latitud<-NULL

dcitri $longitud<-NULL

dcitri$trampa<-NULL

dcitri$Total<-NULL

Ahora hay que transformar el orden de las columnas para apilarlas en una sola, para esto con la funcion melt de reshape2. La siguiente instrucción apila las columnas 1, 2,..., 12 en una sola y con el nombre "MONTH", los valores se apilan en otra columna con el nombre "COUNT".

dcm<-melt(dcitri,variable.name="MONTH",value.name='COUNT')

Ahora procedemos a agrupar los datos por mes:

cats<-group_by(dcm,MONTH)

Calcular Estadísticas Descriptivas e Índices de agregación

Para calcular las estadísticas descriptivas, aplicar summarise con la lista de las estadísticas y su nombre:

dcstats<-summarise(cats,dcmed=mean(COUNT),dcvar=sd(COUNT)^2,ttl=sum(COUNT),ns=length(COUNT))

La instrucción anterior calcula la media dcmed, varianza dcvar y la suma de las observaciones ttl, así como el tamaño de muestra ns de cada mes:

dcstats

   MONTH      dcmed      dcvar ttl ns

1     X1  0.2708333   1.095301  13 48

2     X2  6.1041667 129.584663 293 48

3     X3  3.0000000  54.851064 144 48

4     X4  1.0625000   3.464096  51 48

5     X5  8.5208333 145.701684 409 48

6     X6 15.7500000 810.148936 756 48

7     X7  8.4583333  73.359929 406 48

8     X8  2.8958333   7.286791 139 48

9     X9  2.9166667  10.205674 140 48

10   X10  1.4791667   7.829344  71 48

11   X11  1.4583333   3.913121  70 48

12   X12  1.7916667   7.402482  86 48

Con base en estos resultados, es sencillo calcular algunos índices de agregación y asignar estos valores a nuevas variables en el conjunto de resultados:

relación varianza media: varianza/media:

dcstats$vmId<-(dcstats$dcvar/dcstats$dcmed)

El parámetro k de la binomial negativa, estimado por momentos:

dcstats$k<-dcstats$dcmed^2/(dcstats$dcvar - dcstats$dcmed)

El índice de aglomeración media de Lloyd, calculado por momentos:

dcstats$mcrowd<-(dcstats$dcmed+dcstats$vmId-1)

Los valores tabulados, procesados posteriormente son:

MES MEDIA VAR TOTAL N k vmId MC

FEB 0.3 1.1 13 48 0.1 4.0 3.3

MAR 6.1 129.6 293 48 0.3 21.2 26.3

ABR 3.0 54.9 144 48 0.2 18.3 20.3

MAY 1.1 3.5 51 48 0.5 3.3 3.3

JUN 8.5 145.7 409 48 0.5 17.1 24.6

JUL 15.8 810.1 756 48 0.3 51.4 66.2

AGO 8.5 73.4 406 48 1.1 8.7 16.1

SEP 2.9 7.3 139 48 1.9 2.5 4.4

OCT 2.9 10.2 140 48 1.2 3.5 5.4

NOV 1.5 7.8 71 48 0.3 5.3 5.8

DIC 1.5 3.9 70 48 0.9 2.7 3.1

ENE 1.8 7.4 86 48 0.6 4.1 4.9

El cuadro anterior indica, para k, valores entre 0.1 y 1.9, indicando agregación; la relación varianza media oscila entre 2.5 y 51.5 indicando también agregación. Por otra parte, el índice de aglomeración media (MC) oscila entre 3.1 y 66.2, indicando también poblaciones agregadas.

Gráficas

Una manera eficiente de agrupar las observaciones es mediante las gráficas de caja (boxplot):

dcbx<-ggplot(dcm,aes(x=MONTH,y=COUNT))+geom_boxplot()

La gráfica anterior contiene las etiquetas de los meses como X1, ..., X12, asi que para que sea mas ilustrativo hay que cambiarlas, junto con el título del eje X e Y:

dcbx+scale_x_discrete(labels=c("FEB","MAR","ABR","MAY","JUN", "JUL","AGO","SEP","OCT","NOV","DIC","ENE"))+xlab("MES")+ylab("CONTEOS")+theme_bw()

La gráfica previa indica que las mayores poblaciones ocurren en marzo, junio, julio (los máximos) y agosto. 

Una gráfica alternativa posiblemente menos informativa pero mas utilizada es representar los promedios junto con su intervalo de confianza como una gráfica de linea y puntos y con barras verticales que representen ese intervalo de confianza. Un intervalo de confianza para la media al 95% se calcula como:

media +- EE*t,

Donde media es la media de muestreo, EE es el error estándar de la media (= sqrt(S^2/n)), es decir la raíz cuadrada del cociente de la varianza y el tamaño de muestra, y t es el valor de Student de tablas con un nivel de confiabilidad (1-a), usualmente fijado en 95%, por lo que t = 1.96

Para generar la gráfica es necesario primero añadir una columna con valores numéricos de los meses:

dcstats$M<-c(1,2,3,4,5,6,7,8,9,10,11,12)

Calcular la mitad del intervalo de confianza para cada mes:

dcstats$HCI<-1.96*sqrt((dcstats$dcvar/dcstats$ns))

Ahora se hace la gráfica de linea con puntos:

fci<-ggplot(dcstats,aes(x=M,y=dcmed))+geom_line()+geom_point(size=6)+theme_bw()

Ahora se añaden las lineas verticales del intervalo de confianza: geom_errorbar()

fci<-fci+geom_errorbar(aes(ymax = dcmed + HCI, ymin=dcmed - HCI),width=0.25)

Ahora se pueden modificar las etiquetas de los meses y los ejes:

fci + scale_x_continuous(breaks= c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12), labels= c("FEB", "MAR", "ABR", "MAY", "JUN","JUL", "AGO", "SEP", "OCT", "NOV", "DIC", "ENE")) + xlab("MES") + ylab("CAPTURAS PROMEDIO DE D. CITRI")

La gráfica es la siguiente:

La gráfica de la relación varianza media en función de tiempo es:

jdVM<- ggplot(dcstats, aes(MONTH, vmId)) + geom_bar(stat= 'identity', fill= '#EBEB00', colour='black')

Sin embargo, es necesario cambiar las etiquetas de los meses por sus valores reales y nombrar los ejes apropiadamente:

jdVM + scale_x_discrete(labels=c("FEB", "MAR", "ABR", "MAY", "JUN", "JUL", "AGO", "SEP", "OCT", "NOV", "DIC", "ENE")) + xlab("MES (2011 2012)") + ylab("RELACION VARIANZA MEDIA") + theme_bw()

Lo que genera la siguiente gráfica:

Ley de Potencia de Taylor

La ley de potencia de Taylor (TPL) es un modelo en donde la varianza depende de la media de manera exponencial:

S^2 = am^b, (1)

donde S^2 es la varianza, m es la media y b es un exponente. El valor de b indican el nivel de agregación de las poblaciones:

b= 1 el arreglo es aleatorio,

b < 1 el arreglo es regular y

b > 1 el arreglo es agregado (Taylor, 1961).

El modelo anterior se transforma a uno lineal mediante logaritmos:

log(S^2) + log(a)+ b*log(m), (2)

Con esto se pueden estimar los parámetros log(a) y b mediante regresión lineal simple:

dcTPL<-lm(log(dcstats$dcvar)~ log(dcstats$dcmed))

El análisis de varianza es el siguiente:

anova(dcTPL)

Response: log(dcstats$dcvar)

                   Df Sum Sq Mean Sq F value  Pr(>F)    

log(dcstats$dcmed)  1 36.204  36.204  70.274 7.8e-06 ***

Residuals          10  5.152   0.515                    

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Se observa que sí hay un efecto significativo del logaritmo de la media sobre la varianza (p < 0.01). Los estimadores de los parámetros se obtienen con summary:

summary(TPL)

Call:

lm(formula = log(dcstats$dcvar) ~ log(dcstats$dcmed))

Residuals:

     Min       1Q   Median       3Q      Max 

-1.08522 -0.54493 -0.04926  0.64109  0.90242 

Coefficients:

                   Estimate Std. Error t value Pr(>|t|)    

(Intercept)          1.3290     0.2869   4.632 0.000933 ***

log(dcstats$dcmed)   1.6386     0.1955   8.383  7.8e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7178 on 10 degrees of freedom

Multiple R-squared:  0.8754, Adjusted R-squared:  0.863 

F-statistic: 70.27 on 1 and 10 DF,  p-value: 7.8e-06

El valor de log(a) = 1.329 y el valor de b es 1.6386, con un coeficiente de determinacion (r^2) = 0.87, por lo que se puede decir que el modelo es significativo y con un ajuste razonablemente bueno. El valor de a es exp(1.329) = 3.77, donde exp() es la función exponencial por haberse usado logaritmos naturales. Por tanto, el modelo de Taylor es:

TPL: S^2 = 3.77*(media^1.6386)

La prueba para ver si b es igual a uno se hace con t = (b - 1)/ SE(b) = (1.6386-1)/0.1955 = 3.26 y es altamente significativo:

1-pt(3.26,10)

0.004287068

pt es la función de probabilidad de la distribución de t con 10 grados de libertad, correspondientes a los grados de libertad del error en el modelo de regresión.

De lo anterior se desprende que el modelo de Taylor indica una distribución en agregados de Diaphorina citri, igual que el resto de los índices de agregación.

Para graficar la TPL hay que usar la transformación log(); geom_smooth() produce la linea de regresión con su intervalo de confianza al 95%:

ft<-ggplot(dcstats, aes(x= log(dcmed), y=log(dcvar))) + geom_smooth(method= lm) + geom_point(size= 4)

Para cambiar los nombres de los ejes y el fondo:

ft + xlab("ln(Media)") + ylab("ln(Varianza)") + theme_bw()

Tamaños de Muestra

Una vez que se ajusta el modelo de Taylor, es posible calcular los tamaños de muestra para hacer una estimación de la densidad media. La formula para calcular los tamaños de muestra es:

n= ((t/(mD))^2)(S^2), (3)

donde n es el número de trampas (en este caso) para estimar a la media con cierta precisión D, t es el valor de Student con un nivel de confiabilidad 1-a= 0.95, siendo el valor de t= 1.96, D es la precisión relativa a la media, usualmente de 0.2 a 0.30 como una proporción de la media, m es la media y S^2 es la varianza. Puesto que la TPL es la ecuación (1): S^2 = a m^b, sustituyendo en la ecuación del tamaño de muestra (2) se tiene:

n= ((t/D)^2)(a m^(b-2)), (4)

La ecuación anterior se puede representar como una función en r, el script completo es como sigue:

# Jose Lopez-Collado. Junio 2020

# Requiere ggplot2, reshape2

#====================================================================================

#==================Script para calcular tama?o de muestra   =========================

#====================================================================================

#Genera valores de medias de 1 a 20

mv<-seq(from=1, to=20)

#parametros de la formula de tama?o de muestra, t para confiabilidad, ap y bp son parametros del TPL

tv<-1.96

ap<-3.77

bp<-1.6386

#funcion para calcular los tama?os de muestra, media es el vector de valores de la media, dd es la precision, t es la confiabilidad, ap y bp son los parametros del modelo de Taylor 

f.ssTPL<-function(media, dd, tt, ap,bp)

{

result<-((tt/dd)^2)*(ap*(media^(bp-2)))

}

#calcular tama?os de muestra para D=0.2 y D=0.3

D20<-f.ssTPL(mv,0.2,tv,ap,bp)

D30<-f.ssTPL(mv,0.3,tv,ap,bp)

#convertir a frame:

dfdc<- data.frame(mv, D20, D30)

#Transformar frame a estructura compatible con ggplot, debe estar activa libreria reshape2

mdf<-melt(dfdc, id= c("mv"), variable.name= "PREC", value.name= "N")

#Graficar tama?os de muestra, debe estar activa la libreria ggplot2

ssplot<-ggplot(mdf, aes(mv, N))+ geom_line(aes(colour= PREC), size= 1) + geom_point(aes(shape= PREC, colour=PREC), size=4) + labs(x= "MEDIA", y= "NUMERO DE TRAMPAS") + theme_bw() + theme(text=element_text(size= 18), legend.position= c(0.8, 0.8)) + scale_colour_manual(values= c('black', 'orange'))

# para presentar la imagen en pantalla hay que escribir ssplot en la linea de comando

#==========================Aqui termina el script           =========================

# guardar estas instrucciones en un archivo de texto con la extension R, por 

#ejemplo ssdcD.R

#ejecutarlo desde la linea de comando con la instruccion:

#source("ssdcD.R")

#Para presentar la gr?fica escribir en la linea de comando

#ssplot

#=============Aqui termina el script           ===================

#==================================================================

La gráfica de tamaños de muestra, o número de trampas amarillas necesarias para calcular a la media con una precisión D y una confiabilidad del 95 es:

Se observa que si la precisión se incrementa de D=0.3 a D=0.2 se incrementan los tamaños de muestra, asimismo, si la media estimada tiene valores cercanos a cero, los tamaños de muestra también se incrementan.

Bibliografía

Bliss, C.I. and R.A. Fisher 1953. Fitting the negative binomial distribution to biological data and note on the efficient fitting of the negative binomial. Biometrics 9: 176- 200

Hernandez-Landa, L., J. Lopez-Collado, C.G. García-García, F. Osorio-Acosta, M.E. Nava-Tablada. 2013. Dinamica espacio-temporal de Diaphorina citri Kuwayama (Hemiptera: Psyllidae) en Murraya paniculata (L.) Jack en Cuitlahuac, Veracruz. Acta Zoologica Mexicana. 29(2): 334-345.

Taylor, L. R. 1961. Aggregation, variance and the mean. Nature (London) 189: 732- 735