Come possiamo lappare o sdolcinare per GLM su più specie separatamente?

Question

Sto cercando di eseguire un GLM su più specie diverse nel mio set di dati. Attualmente ho sottoimpostato i miei dati per ogni specie e ho copiato questo codice ed è diventato un bel casino. So che ci deve essere un modo migliore per farlo, (forse con la funzione lapply?) Ma non sono sicuro di come iniziare con quello.

Sto eseguendo il modello sul CPUE (cattura per unità di sforzo) per una specie e utilizzo Anno, Salinità, Scarico e Pioggia come variabili esplicative.

I miei dati sono qui: https://drive.google.com/file/d/1_ylbMoqevvsuucwZn2VMA_KMNaykDItk/view?usp=sharing

Questo è il codice che ho provato. Ottiene il lavoro fatto, ma ho appena copiato questo codice e cambiato la specie ogni volta. Spero di trovare un modo per semplificare questo processo e ripulire un po 'il mio codice.

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)

Produzione:

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

Answer 1

Suggerirei il prossimo approccio creando una funzione per i modelli e quindi utilizzando lapplyun elenco che risulta dall'applicazione split()al dataframe per variabile 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)

In Lmodsci saranno i risultati dei modelli, qui un esempio:

Lmods$`Atlantic Stingray`

Produzione:

 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