我们如何分别在多个物种上申请或申请GLM？

问题描述

我正在尝试对数据集中的多个不同物种运行GLM。目前，我正在为每个物种子集数据并复制此代码，结果变得一团糟。我知道必须有一个更好的方法（也许使用lapply函数？），但是我不确定如何开始。

我正在一个物种的cpuE（每单位工作量捕捞量）上运行该模型，并使用年，盐度，流量和降雨量作为解释变量。

我的数据在这里：https://drive.google.com/file/d/1_ylbMoqevvsuucwZn2VMA_KMNaykDItk/view?usp=sharing

这是我尝试过的代码。它可以完成工作，但是我一直在复制此代码并每次更改种类。我希望找到一种简化此过程并清理代码的方法。

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,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)

输出：

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

解决方法

我建议使用以下方法为模型创建函数，然后在列表上使用lapply，这是通过将split()通过变量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,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)

在Lmods中，将有模型的结果，这里是一个示例：

Lmods$`Atlantic Stingray`

输出：

 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

glm r r