0

I am attempting to model an experiment which was measuring the production of co2 and ch4 at varying temperatures. The experiment has two underlying treatments: substrate with three levels and climate with three levels. I have tried fitting both a gamma glmer with log link and a log transformed lmer. The glmer fails to converge even with optimization, the lmer has no convergence issues is struggling with prediction accuracy at the lower and higher bounds.

#non-converging model

full_model_Gamma <- glmer(co2_accumulation ~ substrate #+ climate 
                      + ch_temp 
                        + treatment
                        + waterbath
                        + ch_temp:substrate
                        + ch_temp:treatment
                        + treatment:substrate 
                        + (1|site) 
                        , family = Gamma(link = 'log')
                        , data = Oxymax#[waterbath == "Old",]
                        , na.action = "na.fail"
                        , control=glmerControl(optimizer="nloptwrap",optCtrl=list(algotithm = "NLOPT_LN_NELDERMEAD")))

#converging model

full_model <- lmer(log(co2_accumulation_s+1) ~ substrate  + ch_temp_s 
               #+ climate  removed from drop1
                        + treatment
                        + waterbath
                        + ch_temp_s:substrate
                        + ch_temp_s:treatment
                        + treatment:substrate # SIR treatments could increase the CO2 rate of calc sites more than period ones due to more microbes
                        #+ substrate:climate
                        #+ (1+treatment|site)
                        + (1|site) 
                        
                        , data = Oxymax#[waterbath == "Old",]
                        , na.action = "na.fail"
                        , control=lmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e5))
                        )

''' plot(full_model) residuasls vs fitted

qqPlot(residuals(full_model, type = "pearson"))

qqplot

plot(residuals(full_model_Gamma, type = "pearson") ~ factor(Oxymax$substrate))

residual distribution across substrates

plot(residuals(full_model_Gamma, type = "pearson") ~ factor(Oxymax$climate))

residual distribution over climate

plot(residuals(full_model_Gamma, type = "pearson") ~ Oxymax$ch_temp) #within this you see potential interaction between ch temp and

residuals over chamber temp

Any help would be greatly appreciated.

Daniel Fishburn
  • 1
  • Why are you using `log(y + 1)` as your transformation? I would also question using the gamma distribution family. – Roland Sep 15 '21 at 07:49
  • Also, your `full_model` is obviously not a "full model" because it excludes higher-order interactions. You seem to expect strong skewness and I wouldn't expect that for CO2 fluxes. – Roland Sep 15 '21 at 07:52
  • @Roland, I forgot to mention that I scaled my response variable so needed to add the plus 1 as it scales form a -1 to +1 range. Also, I understand that it technically isn't a full model but it was the best fitting one. I used drop1 to remove non-significant interactions to improve the model power. – Daniel Fishburn Sep 15 '21 at 08:02
  • I'm sorry but I believe your scaling shouldn't be done (at least the way you do it). I suspect you've just created additional issues for yourself by doing this. "drop1 to remove non-significant interactions" does not "improve the model power", see https://stats.stackexchange.com/a/18245/11849. – Roland Sep 15 '21 at 08:09
  • @Roland I have also ran it without scaling with similar outputs. Thank you for referring me to that post. I have only recently began the journey into modelling so still have much to learn. – Daniel Fishburn Sep 15 '21 at 09:08

0 Answers0