Plotting categorical covariate against occupancy using unmarked package

Question

I am plotting covariates against predicted occupancy models using the unmarked package, three of my covariates are continuous so I have plotted using predict function and ggplot, geom_ribbon. However obviously this doesn't work with categorical/factor variables, I would like to be able to plot as a boxplot the predicted occupancy for the two discrete categories within my factor covariate.

The dataset UMF - is the unmarked frame with site covariates, obs covariates and capture history of the individual species. I have included the code for the null model, continuous covariate model (path_dist) and the point I am at with the categorical covariate (fox_presence). The categorical covariate has two levels: present and abscent and is treated as a factor in the dataset. I have tried to use the same predict function as with the continous and null but changed the type to "response" however this produces an error code.

Is there any way I can model and plot categorical covariates against individual species occupancy in the unmarked package? I have cut out the modelling of the other continuous variables as its just repition but is why the model moves from m2 to m5.

m1 <- occu(formula = ~1
~1,
data = umf) 

m2 <- occu(formula = ~1 # detection formula first 
~path_dist, # occupancy formula second, 
data = umf) 

newDat <- 
cbind(expand.grid(path_dist=seq(min(cov$path_dist),max(cov$path_dist), 
length.out=100)))

newDat<- predict(m2, type="state", newdata = newDat, appendData=TRUE) # 
predict psi (type = "state") and confidence intervals based on our 
model for these road distances 

p1 <- ggplot(newDat, aes(x = path_dist, y = Predicted)) + 
geom_ribbon(aes(ymin = lower, ymax = upper), alpha = 0.5, linetype = 
"dashed") + geom_path(size = 1) + labs(x = "Distance to path", y = 
"Occupancy probability") + theme_classic() + coord_cartesian(ylim = 
c(0,1)) 

#fox_presence model/categorical covariate
cov$Fox_Presence <- factor(cov$Fox_presence) 

m5 <- occu(formula = ~1 
~Fox_presence, data = umf) 

newDat4 <- cbind(expand.grid(Fox_presence=seq(cov$Fox_presence), 
(cov$Fox_presence), length.out=100))
newDat4 <- predict(m5, type="response", newdata = newDat4, 
appendData=TRUE) 

#error: valid types are state, det