¿Cómo se aplica o se aplica la aplicación de GLM en varias especies por separado?

Aug 21 2020

Estoy tratando de ejecutar un GLM en varias especies diferentes en mi conjunto de datos. Actualmente he estado subconfigurando mis datos para cada especie y copiando este código y se ha convertido en un gran lío. Sé que tiene que haber una mejor manera de hacer esto (¿tal vez con la función lapply?), Pero no estoy seguro de cómo comenzar con eso.

Estoy ejecutando el modelo en la CPUE (captura por unidad de esfuerzo) para una especie y usando Año, Salinidad, Descarga y Precipitación como mis variables explicativas.

Mis datos están aquí: https://drive.google.com/file/d/1_ylbMoqevvsuucwZn2VMA_KMNaykDItk/view?usp=sharing

Este es el código que he probado. Hace el trabajo, pero acabo de copiar este código y cambiar la especie cada vez. Espero encontrar una manera de simplificar este proceso y limpiar un poco mi código.

fish_df$pinfishCPUE <- ifelse(fish_df$Commonname == "Pinfish", fish_all$CPUE, 0) #create binomial column fish_df$binom <- ifelse(fish_df$pinfishCPUE > 0, 1,0)


glm.full.bin = glm(binom~Year+Salinity+Discharge +Rainfall,data=fish_df,family=binomial)
glm.base.bin = glm(binom~Year,data=fish_df,family=binomial)

#step to simplify model and get appropriate order
glm.step.bin = step(glm.base.bin,scope=list(upper=glm.full.bin,lower=~Year),direction='forward',
                    trace=1,k=log(nrow(fish_df)))

#final model - may choose to reduce based on deviance and cutoff in above step
glm.final.bin  = glm.step.bin
print(summary(glm.final.bin))

#calculate the LSMeans for the proportion of positive trips
lsm.b.glm = emmeans(glm.final.bin,"Year",data=fish_df)
LSMeansProp = summary(lsm.b.glm)

Salida:

Call:
glm(formula = log.CPUE ~ Month + Salinity + Temperature, family = gaussian, 
    data = fish_B_pos)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8927  -0.7852   0.1038   0.8974   3.5887  

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  2.38530    0.72009   3.313  0.00098 ***
Month        0.10333    0.03433   3.010  0.00272 ** 
Salinity    -0.13530    0.01241 -10.900  < 2e-16 ***
Temperature  0.06901    0.01434   4.811  1.9e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.679401)

    Null deviance: 1286.4  on 603  degrees of freedom
Residual deviance: 1007.6  on 600  degrees of freedom
AIC: 2033.2

Number of Fisher Scoring iterations: 2

Respuestas

1 Duck Aug 21 2020 at 04:56

Sugeriría el siguiente enfoque crear una función para los modelos y luego usar lapplyuna lista que resulta de aplicar split()al marco de datos por variable Commonname:

library(emmeans)
#Load data
fish_df <- read.csv('fish_df.csv',stringsAsFactors = F)
#Code
List <- split(fish_df,fish_df$Commonname) #Function for models mymodelfun <- function(x) { #Create binomial column x$binom <- ifelse(x$pinfishCPUE > 0, 1,0)
  
  
  glm.full.bin = glm(binom~Year+Salinity+Discharge +Rainfall,data=x,family=binomial)
  glm.base.bin = glm(binom~Year,data=x,family=binomial)
  
  #step to simplify model and get appropriate order
  glm.step.bin = step(glm.base.bin,scope=list(upper=glm.full.bin,lower=~Year),direction='forward',
                      trace=1,k=log(nrow(x)))
  
  #final model - may choose to reduce based on deviance and cutoff in above step
  glm.final.bin  = glm.step.bin
  print(summary(glm.final.bin))
  
  #calculate the LSMeans for the proportion of positive trips
  lsm.b.glm = emmeans(glm.final.bin,"Year",data=x)
  LSMeansProp = summary(lsm.b.glm)
  return(LSMeansProp)
}
#Apply function
Lmods <- lapply(List,mymodelfun)

En Lmodshabrá los resultados de los modelos, aquí un ejemplo:

Lmods$`Atlantic Stingray`

Salida:

 Year emmean    SE  df asymp.LCL asymp.UCL
 2009  -22.6 48196 Inf    -94485     94440

Results are given on the logit (not the response) scale. 
Confidence level used: 0.95