## Modelling red fox home range size ## # April 2020 # # Load the appropriate packages library(lme4) library(MuMIn) library(arm) # Import the data and check the structure foxdf<-read.csv("fox_data.csv") str(foxdf) # Create a correlation matrix to identify any predictor variables that may be auto-correlated cor(foxdf[,c("NDVI_mean","NDVI_seas","HFI","Temp_mean","Temp_seas","Precip_mean","Precip_seas")]) # NDVI_seas has a moderately high correlation with Temp_mean and Temp_seas, so we should exclude it from the model # Fit a global model containing the predictor variables and interactions of interest m0<-lmer(HR_mean ~ NDVI_mean * HFI + Temp_mean + Temp_seas + Precip_mean + Precip_seas + (1| Study_ID), data = foxdf, REML=FALSE) plot(m0) # Check the residual plot to see if the model is a good fit for the data (residual plot doesn't look good) # Refit the model with the response variable log-transformed m1<-lmer(log_HR_mean ~ NDVI_mean * HFI + Temp_mean + Temp_seas + Precip_mean + Precip_seas + (1| Study_ID), data = foxdf, REML=FALSE) plot(m1) # Check the residual plot (it looks a lot better) m1_std<-standardize(m1) # Standardise the predictor variables using the arm package # Use the dredge function from MuMIn to fit all possible combinations of predictors from the global model options(na.action = "na.fail") modsel_1<-dredge(m1_std); modsel_1 # Note that the 'boundary (singular) fit' warning occurs because Study_ID accounts for zero variance in some of the models. This is not a problem for model fit or intepretation, though. # Get the parameter estimates and R2 values for the top model. Change the '1' inside get.models() to 2, 3, 4, to get values for the second, third and fourth models (etc). summary(get.models(modsel_1, 1)[[1]]);confint(get.models(modsel_1, 1)[[1]]);r.squaredGLMM(get.models(modsel_1, 1)[[1]]) ################################### # Analyse the HFI components separately # Create a correlation matrix to exclude correlated variables (r > 0.5) from featuring in the same model smat <- abs(cor(foxdf[, c("Built_areas", "Cropland", "Lights", "Nav_waterways", "Pasture", "Human_popn", "Railways","Roads")])) <= 0.5 smat[!lower.tri(smat)] <- NA z.smat<-smat colnames(z.smat) <- paste("z.", colnames(z.smat),sep="") rownames(z.smat) <- paste("z.", rownames(z.smat),sep="") # Fit the model m2<-lmer(log_HR_mean ~ Built_areas + Cropland + Lights + Nav_waterways + Pasture + Human_popn + Railways + Roads + (1| Study_ID), data = foxdf, REML=FALSE) # Note that the 'boundary (singular) fit' warning occurs because Study_ID accounts for zero variance in the global model (see summary(m2)). The same applies when the dredge function is run below. This is not a problem for model fit or intepretation, though. m2_std<-standardize(m2) # Standardize the predictor variables # Use the dredge function to fit all possible combinations of predictors from the global model, excluding models containing highly correlated variables (subset term) modsel2<-dredge(m2_std, subset=z.smat) modsel2_subset<-modsel2[-18,] # Exclude the model that did not converge modsel2_subset<- model.sel(modsel2_subset); modsel2_subset # Get the parameter estimates for the top model summary(get.models(modsel2_subset, 1)[[1]]);confint(get.models(modsel2_subset, 1)[[1]]);r.squaredGLMM(get.models(modsel2_subset, 1)[[1]]) # Compare models that only contain HFI or human population density m_hfp<-standardize(lmer(log_HR_mean ~ HFI + (1| Study_ID), data = foxdf, REML=FALSE)) m_hpd<-standardize(lmer(log_HR_mean ~ Human_popn + (1| Study_ID), data = foxdf, REML=FALSE)) summary(m_hfp);AICc(m_hfp);r.squaredGLMM(m_hfp) summary(m_hpd);AICc(m_hpd);r.squaredGLMM(m_hpd)