I have been following a guide on occupancy modelling provided by Dr. James Paterson:

Occupancy Models in R (Dr. James Paterson)

The global occupancy model I created is as such:

INBU_global_mod <- occu(formula = ~scale(year) + scale(dates) + scale(roads_100m) #detection formula first
                       ~scale(vl_100m) + scale(roads_100m) + scale(build_100m) + year + (1|neigh) + (1|site),
                       data = INBU_stat.sample.unmarkedFrame_cov)

I then created a list of all models, and a ‘newdata’ dataframe:

occu_vl100_newdata <- data.frame(vl_100m = seq(min(INBU_stat.site_cov$vl_100m), 
                                           max(INBU_stat.site_cov$vl_100m), by = 0.23890),
                             vl_50m = mean(INBU_stat.site_cov$vl_50m), # hold other variables constant
                             ndvi_50m = mean(INBU_stat.site_cov$ndvi_50m),
                             ndvi_100m = mean(INBU_stat.site_cov$ndvi_100m),
                             build_50m = mean(INBU_stat.site_cov$build_50m),
                             build_100m = mean(INBU_stat.site_cov$build_100m),
                             roads_50m = mean(INBU_stat.site_cov$roads_50m),
                             roads_100m = mean(INBU_stat.site_cov$roads_100m),
                             imperv_50m = mean(INBU_stat.site_cov$imperv_50m),
                             imperv_100m = mean(INBU_stat.site_cov$imperv_100m),
                             site = as_factor(INBU_stat.site_cov$site),
                             year = as_factor(INBU_stat.site_cov$year),
                             neigh = as_factor(INBU_stat.site_cov$neigh))

Then I attempt to use the following code to create a modavgPred object:

occu_vl100m_pred <- modavgPred(INBU_mods_full,
                           parm.type = "psi", # psi = occupancy
                           newdata = occu_vl100_newdata, re.form = NA)[c("Parameter",
                                                           "Lower.CL",
                                                           "Upper.CL")]

But I receive the following error:

Error in .local(obj, coefficients, ...) :ncol(coefficients) == length(est) is not TRUE

I have a hunch that I need to somehow control for the random effects in my model(s). I already tried to convert them to factors, but the issue persists.

I am hoping to use the following code:

occu_vl100m_pred_df <- data.frame(Predicted = occu_vl100m_pred $mod.avg.pred,
                              lower = occu_vl100m_pred $lower.CL,
                              upper = occu_vl100m_pred $upper.CL,
                              occu_vl100m_pred)

…then ggplot, to create a figure with occupancy probability on the y-axis and the variable vl_100m on the x-axis.

If anyone could provide any sort of guidance, I would be grateful! This has left me stuck for a while haha.

Thank you! 〳 ° ▾ ° 〵