Dealing with quasi-complete separation in General Additive Model?

Question

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: mediafire.com/file/u50jca8p8kwir4j/firedataset.csv/file)

I create a GAM model in R for the species:

g1 <- gam(bird~s(point, bs="re") + 
            bef.aft*ctrl.imp + 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:

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	114
control	after	39
impact	before	74
impact	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.

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.