doubts about emmeans and post hoc comparison in a nesting variable

Question

Good day, I have a doubt about emmeans, Im doing research in two parks, in each park I have a county and in each county, I have three habitats and on each habitat, I have collected beetles (Count response variable).... this is the sampling design for one park (sampling design for CS, the same for CSI): [![sampling design for park CS, the same for CSI][1]][1]

Park<-c("CS","CS","CS","CS","CS","CS","CS","CS","CS","CSI","CSI","CSI","CSI","CSI","CSI","CSI","CSI","CSI")
county<-c("NO","NO","NO","PA","PA","PA","PB","PB","PB","LA","LA","LA","MO","MO","MO","CO","CO","CO")
Hab<-c("Forest","Highland","Grass","Forest","Highland","Grass","Forest","Highland","Grass","Forest","Highland","Grass","Forest","Highland","Grass","Forest","Highland","Grass")
Abu<-c(10,6,4,22,15,3,7,9,6,42,4,3,4,3,2,1,5,13)
data<-as.data.frame(cbind(Park,Sta,Loc,Hab,as.numeric(Abu)))
        

I've made a Negative Binomial GLM and I'm trying to make pairwise comparisons between each factor for each variable. As you all can see this sampling design have the county variable nested within Park, I'm trying to do the comparison between parks only, it gives me that they're pretty different, I think that it has the effect of the nested variable included.

ulm<-glm.nb(Abu~Park+county %in% Park+Hab,data=data)
marginal = emmeans(ulm, ~ Park) ## if I put the nesting=NULL it gives me NA
pairs(marginal,adjust="tukey")
plot(marginal, comparisons = TRUE)

my question is....how can I examine the differences between park without the county effect? with the same model of course. (I know that make separate models is kind of p-hacking thing) Thanks!. ` [1]: https://i.stack.imgur.com/eNtmp.png

Answer 1

First of all, park is nested in county, not the other way around. Look at a design book. It will tell you that A is nested in B if knowing which A also tells you which B. Analogous to your situation, knowing it's Yosemite tells you it's in California.

I suggest just using a model r.h.s. of ~ county + Park + Hab, or perhaps ~(county + Park)*Hab. That will produce fewer dependencies (NA regression coefficients).

It appears that emmeans correctly detected the nesting structure in spite of the way you specified the model. And is designed to show nesting factors when you specify a nested factor.