1
$\begingroup$

I am analysing a data set that is created from walking transects and recording counts for each group size of animals observed. Each transect has 41 repeats, which was approximately 80% zeros. However, within the actual observations, the counts also vary greatly with a large number of small group sizes but some >60 individuals as well.

I want to analyse these data to determine the type of habitat that is preferable to this species by modeling the count data against the environmental covariates that occur at each transect. I have begun with a logistic regression, for which I converted to a logistic vector for all surveys where 1=presence and 0="pseudo absence". The raw count frequencies are as in the screenshot (sorry I didn't know how else to export this without creating a mess of numbers)

enter image description here

I would now like to look at a finer scale of how the preference of habitat varies among the habitats indicated with a higher proportional "selection" from the logistic regression model. Thus I subsetted the data to just the observations (counts>0). Which have the following distribution (x-axis is group size)

enter image description here

I suspected even without zeros these data were overdispersed I read up on this post and after trying the odTest and @BenBolker's own brand I determined they were.

overdisp_fun <- function(model) {
          rdf <- df.residual(model)
          rp <- residuals(model,type="pearson")
          Pearson.chisq <- sum(rp^2)
          prat <- Pearson.chisq/rdf
          pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
          c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}

 overdisp_fun(mNB)
       chisq        ratio          rdf            p 
217168746.02     88351.81      2458.00         0.00 

Where the mNB object was produced from the following model:

mNB<-glmer.nb(size~scale(Acacia.sp)
        + (1|TransectID)
        , glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000))
        , na.action = na.fail
        , weights = w    # a weighting to account for detectability
        , data=imp.NB)

Which always returns the following singularity warning:

boundary (singular) fit: see ?isSingular

Which I am quite sure is bad news! so I was not surprised to see the high level of overdispersion. I have tried a number of combinations to vary the model's complexity. I have also run a glm.nb which was able to converge but I would prefer to include the random effects to account for the multiple samples taken from the same transect.

I would like to ask if anyone out there might have an idea for what I could try next to improve the fit for these data while including the mixed effects and accounting for the overdispersion.

EDIT

After process of elimination, I have found that the culprit causing the singularity issue is the weights argument which indicates a vector with values from 4.350753 to 6.992258. In most cases when the weight is removed, the model converges. Can anyone explain why and possible remedies?

I also have a question relating to this topic here

$\endgroup$
  • 2
    $\begingroup$ link for a similar problem that shows one strategy for overcoming singularity stats.stackexchange.com/questions/30465/… Check the first reply by Michael Chernick. Essentially, you have a situation in which determinant of a square matrix is zero. This is telling you that you have some sort of linear dependency. So, one possible resolution would be to build model by sequentially adding variables, until you hit the problem of singularity. Then you could make a judgement call as to which variable to retain. $\endgroup$ – PsychometStats Oct 24 '19 at 17:02
  • $\begingroup$ @PsychometStats thanks for the suggestion. Unfortunately, I have already tried this. No matter which covariate I use I get a run of the same singularity warning. I have also tried using year and season as random effects to see if that was the problem but it just keeps throwing the same thing. Using vcov for the above model I get values really close to zero for both (5e-5 4e-6)...any ideas? $\endgroup$ – Nebulloyd Oct 24 '19 at 21:02
  • 1
    $\begingroup$ I’d be curious to know what the variance for your random intercepts looks like with and without the weights. $\endgroup$ – Matt Barstead Oct 26 '19 at 23:44
  • 1
    $\begingroup$ So without actually having your data to play with, what seems like is happening is that your weights are accounting for 100% of your random effect variance, and so after weighting, you have no random intercept variance left for the model to chew on. $\endgroup$ – Matt Barstead Oct 27 '19 at 23:03
  • 1
    $\begingroup$ What exactly goes into your weight value? $\endgroup$ – Matt Barstead Oct 27 '19 at 23:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.