I am modelling influence of fire on occurrence of certain bird species (count response variable) in Before/after control/impact experiment design. I've got intact and burned sites and count data both from before and after the burning. I want to see if the fire has an effect on species occurrence. I've got three random variables: site (18 lvls), year of count (3 lvls) and point count spot within site (10 counts per 18 sites = 180 lvls). I use family=nb based on my data distribution (histogram, slight overdispersion).

dataset fragment as csv file

I create a GAM model in R for the species:

mydata <- read.table("clipboard",header=TRUE,sep="\t",na.strings="NA",dec=".",strip.white=TRUE)
mydata$point.count.id <- as.factor(mydata$point.count.id)
mydata$site<- as.factor(mydata$site)
mydata$year <- as.factor(mydata$year)
mydata$bird.numbers  <- as.numeric(mydata$bird.numbers)
mydata$before.after <- as.factor(mydata$before.after)
mydata$before.after <- relevel(mydata$before.after, "before")
mydata$ctrl.impact <- as.factor(mydata$ctrl.impact)
g1 <- gam(bird.numbers~s(point.count.id, bs="re") + 
            before.after*ctrl.impact + s(site, bs='re') + s(year, bs='re'),
          data=mydata, family=nb, method="REML")

When I run the model (which seem to work well for other species), I get enormous standard error for the interaction parameter. This is the exact result:

enter image description here

I investigated my data and I think I've found the root of the problem. There was not a single observation of the species on sites struck by fire after the fire (what I believe is a case of "quasi-complete separation"). The summary table for my data looks like this:

impact period sum
control before 117
control after 39
fire before 74
fire after 0

When I experimentally changed one point count observation to value "1" to make impact/after total sum="1" instead of "0", the model results changed drastically and the fire turned out to have a visible impact on the species, which is the expected result.

I've looked up for some advice how to deal with the problem, but mostly found solutions relevant for situations where this perfect prediction is undesired. In my case it's different - this is my most important result.

Is there any way to deal with this problem? If it requires switching to different model than GAM, I'm fine with that (my explanatory variables are super simple). However, GAM-based solution would be the most desired one (for consistency in future publication).

Any advice would be highly appreciated.


Tried fitting the data (response as binomial) into model fitted by Penalized ML (following advices here: How to deal with perfect separation in logistic regression?):

gfirth <- logistf(bird.numbers > 0 ~ before.after*ctrl.impact, data=mydata)

This actually worked quite great, but unfortunately it doesn't take into consideration the random effect. It would be amazing to have a solution which includes random effects.

Other methods I've tried so far:


confint(glmer.nb(bird.numbers~before.after*ctrl.impact + (1|point.count.id) +
 (1|year) + (1|site), data=mydata)) 

CI still enormous


confint(glm(bird.numbers~before.after*ctrl.impact, data=mydata,
         family = negative.binomial(theta = 48417.29, link = "log")))

CI still enormous. [Theta value derived from previous glmer.nb model.]


g1 <- gam(bird.numbers~s(point.count.id, bs="re") + 
            before.after*ctrl.impact + s(site, bs='re') + s(year, bs='re'),
          data=mydata, family=nb, method="REML")
gam.mh(g1, ns=1000,burn = 100, t.df=40, rw.scale=.5,thin=1)

Getting this error:

"in if (u[i] < exp(lfp[j] - lfp[i] + lpl1 - lpl0)) { : missing value where TRUE/FALSE needed"

Tried experimenting with parameter values (I do not understand them quite well), but nothing worked. Changing the species doesn't result in an error.


g2 <- gam(bird.numbers~s(point.count.id, bs="re") + 
            before.after*ctrl.impact + s(site, bs='re') + s(year, bs='re'),
          data=mydata, family=nb, method="REML", fit=FALSE)
ginla(g2, A=NULL, nk=16, nb=100, J=1, interactive=FALSE, int=0, approx=0)

Getting this error:

"in na.fail.default(data.frame(x = as.numeric(obj1), y = as.numeric(obj2))) : missing values in object In addition: Warning message: In min(dens0, na.rm = TRUE) : no non-missing arguments to min; returning Inf"

Again tried experimenting with parameter values (again, I do not understand what they mean too well), but nothing worked. Changing the species doesn't result in an error.

Right now trying to figure out how bootstrapping methods works, but couldn't work up to a code yet - it's quite complicated to me so far.

Perhaps I messed something up with ginla or gam.mh parameters? Could anyone suggest how to estimate parameter using bootstrapping methods? Maybe anyone else has some other suggestions to try out, which would include random effects?

  • $\begingroup$ Similar: stats.stackexchange.com/questions/400893/… $\endgroup$ Commented Jan 14, 2021 at 17:59
  • $\begingroup$ Some other things to clean up in the post: 1) R package used should be loaded, so the code should iinclude library(mgcv) and maybe others ... 2) variable names in code should match names in data file, and (re)coding such as declaring as factors should be done in the included code $\endgroup$ Commented Jan 24, 2021 at 19:37

1 Answer 1


Your diagnosis is probably correct, and the big standard errors is because the loglikelihood function in this case is far from quadratic. Confidence intervals based on likelihood profiling might be better, as discussed at Why is there a difference between manually calculating a logistic regression 95% confidence interval, and using the confint() function in R? and elsewhere (using the confint function in R.) Unfortunately, there does not seem to be a confint method for gam objects from mgcv. Maybe some other R package can be ised for profiling in that case (if you could post your dataset we could try.)

Some other ideas: Try the gam.mh function for posterior sampling, for a Bayesian analysis based on your estimated model object. Or parametric bootstrap ... or the ginla function.

Edit Some more comments: Trying to run your analysis, but tere are multiple problems. Please update code and data file so variable names are the same, and post complete R code including loading of used R packages, any transformations/name changes, so your code can be run directly from the imported .csv file. That cannot be done now.

Also, check that point and year are coded directly, their estimated variances are close to zero, and the deviance of the two models with/without are almost the same.

  • 1
    $\begingroup$ Thank you kindly for your answer. I am submitting the dataset fragment in case someone wants to help/practice. Unfortunately, the dataset contains 540 rows and i wasn't sure how to create a short representative one, so i submit original dataset fragment as a csv on mediafire (see my original post). I have an obligatory break until 21st jan, so before i ask kindly for any further help, i shall first try tackling the problem myself following your suggestions. I will update this Q and let know of my results. Thank you. Best, $\endgroup$ Commented Jan 18, 2021 at 12:46
  • $\begingroup$ EDIT: I've updated original post according to your suggestions and included all data transformation and packages - now it should run smoothly. Sorry for the mess. $\endgroup$ Commented Jan 25, 2021 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.