Random (Selection) Skewers.

Volví de mis vacaciones con algo pesado…

Dos poblaciones sujetas a un escenario idéntico de selección natural pueden diferir en la forma que responden porque difieren en G, la matriz de varianza-coavarianza genética entre los caracteres. La "forma" de G (en general se representa como una pelotita, más o menos deforme, más alargada en la dirección donde hay más variabilidad) puede acelerar o retardar el cambio evolutivo.

Invirtiendo el razonamiento, podemos comparar las Gs de dos poblaciones. Esta es la idea básica de los random (o selection) skewers propuestos por Cheverud (1-aquí una buena descripción). Como no siempre tenemos datos de la selección, la simulamos creando miles de vectores de selección, aplicando cada uno de ellos (skewer se traduce como "brocheta") a ambas matrices G, obteniendo la dos respuestas a la selección y evaluando si estas respuestas difieren en su magnitud y su dirección.

Supongamos que tenemos dos matrices G (G1 y G2), con el mismo número de caracteres N

N<-ncol(G1)

Creamos un vector de largo N extraido de una distribución uniforme

RS1<-runif(n=N)

RS1

Pero necesitamos que contega valores positivos y negativos con el 50% de probabilidad para cada uno

sign<-c(-1,1)

RS2<-vector(length=N)

for (i in 1:N) RS2[i]<-sample(sign, size=1)*RS1[i]

RS2

Y ahora necesitamos que tenga una longitud estándar, específicamente que la suma de sus elementos al cuadrado sume 1.

RS3<-RS2/sqrt((t(RS2)%*%RS2))

RS3

El resto ya es sencillo. Obtenemos las dos respuestas a la selección (dz)

dz1<-G1%*%RS3

dz2<-G2%*%RS3

Comparamos su magnitud

L<-sum(dz1^2)/sum(dz2^2)

L

Comparamos su correlación (el coseno de su ángulo en el hiperespacio fenotípico, ja!)

A<-(t(dz1)%*%dz2)/ sqrt((t(dz1)%*%dz1)*(t(dz2)%*%dz2))

A

Ahora todo en una sola función…

selection.skewers<-function(G1, G2){

N<-ncol(G1)

RS1<-runif(n=N)

sign<-c(-1,1)

RS2<-vector(length=N)

for (i in 1:N) RS2[i]<-sample(sign, size=1)*RS1[i]

RS3<-RS2/sqrt((t(RS2)%*%RS2))

dz1<-G1%*%RS3

dz2<-G2%*%RS3

L<-sum(dz1^2)/sum(dz2^2)

A<-(t(dz1)%*%dz2)/ sqrt((t(dz1)%*%dz1)*(t(dz2)%*%dz2))

res<-c(L, A)

res}

1- Cheverud J.M. y Marroig G. (2007). Comparing covariance matrices: random skewers method compared to the common principal components model. Genetics and Molecular Biology, 30: 461–469.

Y entonces resulta que la selección correlacional era significativa…

Diantres, y ahora debería hacer una superficie de selección para ver que pasa…

Bueno, hay varias opciones, pero como sabemos que las superficies paramétricas (esas que resultan de usar los parámetros del modelo de selección de Lande y Arnold), pueden dar una idea engañosa de la realidad y sólo son capaces de representar un rango limitado de formas (no lo digo yo…) vamos a usar una superficie no paramétrica, cuyo origen se rastrea a Dolph Schluter y sus splines cúbicos (1,2).

Hace un tiempo, leyendo el paper de Mcglothlin et al. (3) encontré cómo hacerlo en R con el paquete fields. Desde entonces lo he aplicado, publicado, y la gente parece no quejarse al respecto.

Necesitan dos medidas fenotípicas y una medida de fitness. la rutina es en extremo sencilla.

#cargamos el paquete

library(fields)

#no deben existir datos faltantes.

datos<- na.omit(datos)

#unimos las dos variables fenotípicas

sup<-cbind(datos$X1, datos$X2)

#ajustamos el modelo

#el primer término son las variables predictoras,

#el segundo es el fitness

fit<-Tps(sup, datos$fitness)

Y ya. La explicación de lo que el programa hizo, es un poquito oscura (tomen el curso de Métodos en Ecología Evolutiva que damos con Mariano Ordano y Andrea Cocucci!! pero no prometo que la entiendan de todas formas). Se ha calculado una aproximación a una superficie no paramétrica, que es lo que necesitamos para visualizar la selección correlacionada. Esto último se puede hacer con image y persp. 

image(predict.surface(fit))

persp(predict.surface(fit))

A lo cual, por supuesto, les podemos manipular los parámetros gráficos para que se vean mejor.

image(predict.surface(fit), col = gray(seq(0, 1, by=0.02)),

xlab="nectario", ylab="número de flores")

contour(predict.surface(fit), add=TRUE, col="red")

points(x=datos$nectario, y=datos$flores, pch=20)

persp(predict.surface(fit), border="slateblue4", col= "aquamarine2", xlab= "nectario", ylab= "número de flores", zlab= "polinarios", theta= -45, phi= 20, r=5, d=1, ticktype= "detailed", shade= 0.1 )

1-Schluter D. 1988. Estimating the form of natural selection on a quantitative trait. Evolution 42: 849-861.

2-Schluter D, Nychka D. 1994. Exploring fitness surfaces. The American Naturalist 143: 597-616.

3-Mcglothlin JW, Parker PG, Nolan VJ, Ketterson ED. 2005. Correlational selection leads to genetic integration of body size and an attractive plumage trait in dark-eyed juncos. Evolution 59: 658-671.

Volviendo a las bases

Algo sencillo, fácil y que causa dolores de cabeza: ¿qué significan los coeficientes de una tabla summary de un modelo lineal con interacción?

Supongamos que tenemos dos variables independientes, una es continua (X) y la otra es categórica (Z) con tres niveles (a, b, c). Si queremos probar la interacción nuestro modelo será 

fit<-lm(Y~X*Z)

y un ejemplo de tabla de salida podría ser...

summary(fit)

---

Coefficients:

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

(Intercept)   0.1078     0.3136   0.344 0.732958   

X             0.9714     0.2561   3.793 0.000505 ***

Zb            2.6422     0.4435   5.957 5.94e-07 ***

Zc            0.2278     0.4435   0.514 0.610429   

X:Zb         -1.5674     0.3622  -4.328 0.000101 ***

X:Zc         -0.7103     0.3622  -1.961 0.057013 . 

---

Residual standard error: 0.7013 on 39 degrees of freedom

Multiple R-squared: 0.6045,     Adjusted R-squared: 0.5538

F-statistic: 11.92 on 5 and 39 DF,  p-value: 4.884e-07

Antes de seguir, extraigamos los coeficientes.

B<-fit$coeff

Esta forma de presentar los resultados se debe a la codificación como dos dummies de la variable categórica, que toman valores de (0;0) para el nivel a (tomado como referencia alfabéticamente), (1;0) para el nivel b y (0;1) para el nivel c. Así, la fórmula completa arrojada por summary hace referencia a

Y = B[1] + B[2]X + B[3]dummy1 + B[4]dummy2 + B[5]Xdummy1 + B[6]Xdummy2.

La tabla de summary tiene toda la información necesaria para dibujar 3 rectas, una por cada nivel del factor. 

Para el nivel a

beta1=0.1078                    beta2= 0.9714

B1a<- B[1];              B2a<- B[2]

Para el nivel b

beta1 =0.1078 + 2.6422      beta2= 0.9714 + (-1.5674)

B1b<- B[1]+B[3];          B2b<- B[2]+B[5]

Para el nivel c

beta1 =0.1078 + 0.2278      beta2= 0.9714 + (-0.7103)

B1c<- B[1]+B[4];          B2c<- B[2]+B[6]

 

Podemos verlo en un gráfico

plot(Y~X, type="n")

abline(B1a,B2a)

abline(B1b,B2b, col="red")

abline(B1c,B2c, col="blue")

 

Finalmente, una tabla anova nos indicará si el término de interacción es globalmente significativo.

anova(fit)

Analysis of Variance Table

Response: Y

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

X          1  1.0129  1.0129  2.0593 0.159250   

Z          2 19.0662  9.5331 19.3825 1.43e-06 ***

X:Z        2  9.2397  4.6199  9.3930 0.000468 ***

Residuals 39 19.1817  0.4918                    

---

Ver como abeja

Una amiga quería clasificar las flores de su planta preferida (véase infierno en Wikipedia) de acuerdo a cómo son vistas por las abejas. Las abejas tienen visión tricromática en verde-azul-ultravioleta, así que para reconstruir cómo ven se necesita un espectrómetro sensible al UV. De allí obtendremos datos con dos columnas: una especifica la longitud de onda en nanómetros y la restante la reflectancia del objeto (la proporción de luz emitida que es reflejada).

En este post reconstruyo en R obtener las coordenadas para utilizar el hexágono de color de Apis mellifera, uno de los modelos más usados de la visión en estos insectos, a partir de lo expuesto en Chittka y Kevan (2005). Este capítulo puede descargarse de la página personal del autor.

Para llevar a cabo esta rutina necesitarán guardar en un archivo los valores de la tabla 4.6 de ese capítulo. En mi caso lo hice utilizando como nombres de columna nm, uv, blue, green, D65 y Leaf. Además necesitan las mediciones del espectrómetro, entre los 300 y 700 nm.

El principal problema a resolver es que los valores tabulados de la sensibilidad de los receptores de color y la irradiancia de la luz del día se encuentran en intervalos de 5 nm, mientras que los intervalos de los espectros de las flores dependen del aparato que usemos para tomarlos (en el caso de mi amiga, 0.22 nm). Por lo tanto hay que realizar una interpolación para calcular los valores faltantes. Esto se puede hacer con approxfun del paquete stats. Los demás pasos son sencillos y se encuentran detallados en la bibliografía.

#FUNCION PARA CALCULAR COORDENADAS EN EL HEXÁGONO

#Los parámetros de esta función son AP, NM y REF

#AP es el marco de datos con la sensibilidad de los 3 
#receptores de la abeja, la irradiancia de la luz del día y 
#la reflectancia de una hoja verde típica

#NM es un vector con los nm tomados con el espectrómetro

#REF es un vector con es la reflectancia medida entre 0 y 1

hexagon.coor<-function(NM, REF, AP){

UVfun<-approxfun(AP$nm, AP$uv)

BLfun<-approxfun(AP$nm, AP$blue)

GRfun<-approxfun(AP$nm, AP$green)

D65fun<-approxfun(AP$nm, AP$D65)

LEAfun<-approxfun(AP$nm, AP$Leaf)

Ru<-1/sum(UVfun(NM)*D65fun(NM)*LEAfun(NM), na.rm=T)

Rb<-1/sum(BLfun(NM)*D65fun(NM)*LEAfun(NM), na.rm=T)

Rg<-1/sum(GRfun(NM)*D65fun(NM)*LEAfun(NM), na.rm=T)

Pu<-Ru*sum(UVfun(NM)*D65fun(NM)*REF, na.rm=T)

Pb<-Rb*sum(BLfun(NM)*D65fun(NM)*REF, na.rm=T)

Pg<-Rg*sum(GRfun(NM)*D65fun(NM)*REF, na.rm=T)

Eu<-Pu/(Pu+1)

Eb<-Pb/(Pb+1)

Eg<-Pg/(Pg+1)

Xh<-0.8667*(Eg-Eu)

Yh<-Eb-0.5*(Eg+Eu)

res<-cbind(Xh, Yh)

res}

Suponiendo que queremos aplicar esta función a un gran número de columnas, necesitaremos usar apply, especificando los parámetros fijos: NM y AP deben contener los datos ya especificados.

DAT<- …  #mediciones de reflectancia

res<-t(apply(X=DAT, MARGIN=2, FUN=hexagon.coor,  NM=…, AP=…))

res

 

#FUNCIÓN PARA GRAFICAR

hexagon.plot<-function(x, box="n", axes=F,...){

Xcor<-c(-0.8667, -0.866, 0, 0.8667, 0.8667, 0)

Ycor<-c(-0.5, 0.5, 1, 0.5, -0.5, -1)

plot(Xcor, Ycor, pch="", xlab="", ylab="", bty=box, axes=axes,

xlim=c(-1,1), ylim=c(-1.1,1.1))

polygon(Xcor, Ycor)

text(x=c(-0.95, 0, 0.95), y=c(-0.5, 1.1, -0.5),

labels=c("UV", "B", "G"))

points(x[,1],x[,2], ...)}

hexagon.plot(res)

El resultado final es un gráfico como el que se muestra abajo, donde la distancia indica el parecido entre dos colores, desde la percepción de Apis mellifera.

Referencia:

Chittka, L. & Kevan, P.G. (2005) Flower colour as advertisement. In Dafni, A., Kevan, P.G., Husband, B.C. (eds.) Practical Pollination Biology. EnviroquestLtd., Cambridge, ON, Canada, pp. 157 – 196.