Comprobación de la realidad mediante el MLG

Question

Necesito un baño de realidad. Tengo un conjunto de datos, donde sé cuántas mariposas individuales de dos especies co-ocurren en un prado (no siempre, sin embargo). Tengo variables adicionales, por ejemplo prado húmedo/seco, intensamente cultivado/no cultivado, porcentaje de superficie alrededor del prado cubierta por prados húmedos o secos...

Es la cabeza del conjunto de datos. 50 filas en total, 25 por especie. Observe que todas las columnas son idénticas excepto el recuento y la especie, lo que indica que proceden del mismo lugar de muestreo.

> head(dej)
  count     type1 type2 perc.for.100m perc.dry.100m perc.wet.100m species
1     1 intensive   dry        13.836        22.724         0.000   reali
2     3 extensive   wet         6.877         1.613        52.213   reali
3     4 intensive   wet        22.770         0.537        44.901   reali
4     6 intensive   dry        17.346        42.322         6.359   reali
5     1 extensive   wet        34.854         9.091        11.950   reali
6     2 extensive   dry        50.387        19.245         0.000   reali
...
26     0 intensive   dry        13.836        22.724         0.000 sinapis
27     0 extensive   wet         6.877         1.613        52.213 sinapis
28     0 intensive   wet        22.770         0.537        44.901 sinapis
29     0 intensive   dry        17.346        42.322         6.359 sinapis
30     1 extensive   wet        34.854         9.091        11.950 sinapis
31     1 extensive   dry        50.387        19.245         0.000 sinapis
...

Me interesa saber si alguna de estas variables influye en las especies y sus respectivos recuentos.

Y este es el resultado del modelo "completo".

glm(formula = count ~ type1 + type2 + perc.for.100m + perc.dry.100m + 
    perc.wet.100m + species, family = poisson, data = dej)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.8458  -1.1414  -0.4546   0.8297   2.2145  

Coefficients:
               Estimate Std. Error z value Pr(>|z|)   
(Intercept)    0.028129   0.523509   0.054  0.95715   
type1intensive 0.196699   0.191960   1.025  0.30551   
type2wet       0.071841   0.334286   0.215  0.82984   
perc.for.100m  0.003741   0.008277   0.452  0.65130   
perc.dry.100m  0.010952   0.010750   1.019  0.30829   
perc.wet.100m  0.007467   0.011596   0.644  0.51960   
speciessinapis 0.597837   0.187689   3.185  0.00145 **

¿Le parece éste el enfoque correcto?

Información adicional

Como nota al margen, basándome en mi exploración de los datos, esperaría que el recuento dependiera (al menos) de la variable type2, por desgracia no es lo que obtuve.

Utilizando la "lógica inversa", he probado si se pueden predecir las especies utilizando mis datos, lo que presumiblemente confirma los resultados anteriores.

Call:
glm(formula = species ~ type1 + type2 + perc.for.100m + perc.dry.100m + 
    perc.wet.100m + count, family = binomial, data = dej)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.6322  -1.0136  -0.1568   1.0592   1.6407  

Coefficients:
                Estimate Std. Error z value Pr(>|z|)  
(Intercept)    -0.351192   1.658052  -0.212   0.8323  
type1intensive -0.170583   0.651611  -0.262   0.7935  
type2wet       -0.107377   1.078726  -0.100   0.9207  
perc.for.100m  -0.002806   0.026807  -0.105   0.9166  
perc.dry.100m  -0.010227   0.036982  -0.277   0.7821  
perc.wet.100m  -0.006486   0.038071  -0.170   0.8647  
count           0.345036   0.153811   2.243   0.0249 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 69.315  on 49  degrees of freedom
Residual deviance: 63.431  on 43  degrees of freedom
AIC: 77.431

EDITAR 1

Aniko se ha dado cuenta de que puede haber una interacción entre el tipo2 y la especie. ¡Ya lo creo!

Call:
glm(formula = count ~ type1 + type2 * species + perc.for.100m + 
    perc.dry.100m + perc.wet.100m, family = poisson, data = dej)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.0859  -1.1350  -0.1947   0.7109   2.7470  

Coefficients:
                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -0.357165   0.559987  -0.638  0.52360    
type1intensive           0.196699   0.191960   1.025  0.30551    
type2wet                 0.704769   0.429087   1.642  0.10049    
speciessinapis           1.145132   0.306847   3.732  0.00019 ***
perc.for.100m            0.003741   0.008277   0.452  0.65130    
perc.dry.100m            0.010952   0.010750   1.019  0.30829    
perc.wet.100m            0.007467   0.011596   0.644  0.51960    
type2wet:speciessinapis -0.962811   0.394038  -2.443  0.01455 *  

EDITAR 2

Después de eliminar los términos no significativos (suponiendo que he encontrado el máximo global por el bien del dragado de datos) la historia da otro giro en la dirección correcta.

Call:
glm(formula = count ~ type2 * species, family = poisson, data = dej)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7080  -1.1617  -0.1582   0.6979   3.1599  

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)               0.1542     0.2673   0.577  0.56408    
type2wet                  0.6821     0.3237   2.107  0.03508 *  
speciessinapis            1.1451     0.3068   3.732  0.00019 ***
type2wet:speciessinapis  -0.9628     0.3940  -2.443  0.01455 *

Answer 1

Siguiendo la respuesta de Nick Sabbe, he aquí la solución GLMM más sencilla que se me ha ocurrido:

dej$location <- factor(rep(1:25,2))
library(lme4)
glmer(count ~ type1 + type2*species + 
   perc.for.100m + perc.dry.100m + perc.wet.100m + 
   (1|location), family = poisson, data = dej)

También sería buena idea comprobar si hay sobredispersión.

Para tu foto, yo intentaría

library(ggplot2)
ggplot(dej,aes(x=type2,y=count))+stat_sum(aes(size=..n..))+
   facet_grid(.~species)

principalmente por la ventaja de stat_sum que mostrará fácilmente dónde hay un exceso de trazado (para simplificar, puede probar con el jittering).

Answer 2

Usted mismo indica que sus mediciones no son independientes (mide la abundancia de ambas especies en los mismos lugares). Por ello, debe corregir las mediciones repetidas.

Prueba lmer del paquete lme4.