I am an early career ecologist/biologist and I have been working on a Generalized Linear Mixed Model in R to determine influential variables to wildlife behavior. I have been learning as I go and I believe this statistical method has a lot of potential for my research but I would like another set of more experienced eyes (or multiple) to make sure I don't have glaring issues that could undermine the validity of my model, particularly on the stats/coding side. Chat GPT has been useful in the learning process but I don't trust it enough for the final product lol!
My research question is about what variables are most influential for determining bioacoustic (animal sounds) activity for different groups of animals (avian, anuran, etc) in an understudied region. I have various sites, each of which had 3 plots in specific habitat types (2 forest types, 1 unforested per site). At each plot, I collected continuous audio recordings which I processed with a deep learning model to detect and sort bioacoustic sounds into 7 animal groups. I am using the count of the detections for each animal group in a model to help identify whether habitat type, size of trees, and/or weather impact the prevalence of detections for each group and which, if any of the variables, seem to be the best predictor of bioacoustic detections. I have the hour of day each detection occurred. For each plot, I collected size and density of trees in the two forest habitats and hourly temperature and humidity data. I've attached my code with a more detailed explanation of my process below.
# Removing na values. Want to remove na but not entire variables
# For structural (height, canopy, diameter, density),
# replace na with 0
Merged_Data$Height.Average[is.na(Merged_Data$Height.Average)] <- 0
Merged_Data$Canopy.Average[is.na(Merged_Data$Canopy.Average)] <- 0
Merged_Data$Average.Density[is.na(Merged_Data$Average.Density)] <- 0
Merged_Data$Diameter.Average[is.na(Merged_Data$Diameter.Average)] <- 0
# For temp/RH replace with mean
Merged_Data$Average.Humidity[is.na(Merged_Data$Average.Humidity)] <-
mean(Merged_Data$Average.Humidity, na.rm=TRUE)
Merged_Data$Average.Temperature[
is.na(Merged_Data$Average.Temperature)] <-
mean(Merged_Data$Average.Temperature, na.rm=TRUE)
# Group detections by Sound Class
Merged_Data %>%
group_by(Sound.Class)
# Factor variables as factors
# Merged_Data$Sound.Class <- as.factor(Merged_Data$Sound.Class)
# Not necessary now that grouped by sound class?
Merged_Data$Site <- as.factor(Merged_Data$Site)
Merged_Data$Habitat <- as.factor(Merged_Data$Habitat)
Merged_Data$Hour <- as.factor(Merged_Data$Hour)
# Scale numeric variables using mutate to equally weight
# each variable in model
Merged_Data <- Merged_Data %>%
mutate(
Height.Average_Scaled = scale(Height.Average),
Diameter.Average_Scaled = scale(Diameter.Average),
Canopy.Average_Scaled = scale(Canopy.Average),
Average.Density_Scaled = scale(Average.Density),
Average.Humidity_Scaled = scale(Average.Humidity),
Average.Temperature_Scaled = scale(Average.Temperature)
)
# Ensure numeric variables are numeric
str(Merged_Data)
Merged_Data$Slope.Average <- as.numeric(Merged_Data$Slope.Average)
Merged_Data$Height.Average_Scaled <-
as.numeric(Merged_Data$Height.Average)
Merged_Data$Canopy.Average_Scaled <-
as.numeric(Merged_Data$Canopy.Average)
Merged_Data$Diameter.Average_Scaled <-
as.numeric(Merged_Data$Diameter.Average)
Merged_Data$Average.Density_Scaled <-
as.numeric(Merged_Data$Average.Density)
Merged_Data$Average.Humidity_Scaled <-
as.numeric(Merged_Data$Average.Humidity)
Merged_Data$Average.Temperature_Scaled <-
as.numeric(Merged_Data$Average.Temperature)
# combine related variables into single predictor (PCA),
# Sound.Class as fixed effect
pca <- prcomp(Merged_Data[, c("Height.Average_Scaled",
"Diameter.Average_Scaled", "Canopy.Average_Scaled")],
scale. = TRUE)
Merged_Data$PC1 <- pca$x[, 1]
# transform detection data
Merged_Data$Log_Detection <- log(Merged_Data$Detection + 1)
# create GLMM with data grouped by sound class with interaction
# btwn temp and RH
GLMMLogPCA <- lmer(Log_Detection ~ Habitat + PC1 +
Average.Density_Scaled +
Average.Temperature_Scaled * Average.Humidity_Scaled + (1|Site) +
(1|Hour), data = Merged_Data,
control=lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 1e5)))
Currently, I have the data grouped by animal group to make sure the number of detections line up with the corresponding animal group. I'm not sure if there are any potential downsides to this grouping method. I have site and time of day as random effects, and the rest are fixed. Because of potential multicollinearity issues, I combined my structural variables into a single principal component analysis predictor. Temperature and humidity are obviously also related but I used an interaction between the two to handle that because I wasn't getting high VIF values for the two.
I replaced na values with 0 for structural data because those were wetland habitats with no trees. I replaced the mean for missing temp and humidity because some data were lost there on account of technical difficulties. I scaled the variables since they were all on different scales and used a log transformation on the detections for normalization purposes.
I also checked for singularity and assumption validation. No singularity but I have encountered some plots (attached below) that I'm not sure how to interpret for the assumptions.
Q-Q plot looks fairly normal, and was improved by the log transformation. Heteroscedasticity plot has a strange diagonal limit. Wondering if this is a result of the 0s in my data and the transformations but haven't been able to find a similar example anywhere to figure out if it's something to worry about or not. I've never used an influence plot before so I'm not sure how worried I should be about the irregular appearance of this one. My data has lots of outliers but they are important so I can't throw them out.
My plan after finalizing the model was to run some Tukey tests on the most influential variables to identify direction and magnitude of their impact.
log(y + c)
transformation with arbitraryc
also indicates that you should probably use a different type of regression model (like maybe a GLMM with different link function and distribution family, possibly even a model accounting for zero-inflation). You need to describe your dependent variable more clearly and in more detail. Is this perhaps a count? $\endgroup$NA
with zero. Tree height isn't zero. If you have missing values, the first question is why they are missing. Only then you can decide how to handle them. That's another separate question. $\endgroup$