¿Cómo se aplica o se aplica la aplicación de GLM en varias especies por separado?
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
Sugeriría el siguiente enfoque crear una función para los modelos y luego usar lapply
una 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 Lmods
habrá 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