3
$\begingroup$

I was trying to work out this on my own but I find myself overwhelmed.

I am looking at the amount of moisture present in forest fuels (FMC) following forest thinning. Forest fuels are classified by their size into five distinct size classes (FuelType). I monitored changes in FMC over the course of one dry season by monthly sampling (SamplingMonth) of each FuelType in 13 different forest stands (UnitID). Forest stands include 10 thinned forest, and 3 unthinned forests (Treatment). Forest stands vary based on their age, and time since thinning took place.

I would like to determine if there are differences in FMC between thinned and unthinned forests. Due to the fact that FuelType respond to seasonal changes differently, it is necessary to consider them separately. Using lme4 package I built the following models:

model.one = lmer(log(FMC) ~ Treatment + FuelType + (SamplingMonth | UnitID) , data=moisture)
model.two = lmer(log(FMC) ~ Treatment + FuelType + (1+SamplingMonth | UnitID) , data=moisture)

I was able to run them both, which is a success in itself. But I am not convinced that these models are really telling me what I wanted to know. For some reason I anticipated a random effect (SamplingMonth | FuelType) in there somewhere but the model will not budge.

Upon request I am attaching a link to my data:
http://radekonline.com/download/ForAnalysis/MoistureDataBasic.dat
And the R code to retrieve it:

moisture <- read.csv("http://radekonline.com/download/ForAnalysis/MoistureDataBasic.dat")

Any help and opinion would be greatly appreciated!

Thank you.

$\endgroup$
2
$\begingroup$

Thanks for the question!

First off, the two models you show are identical - (SamplingMonth | UnitID) is implicitly (1 + SamplingMonth | UnitID). If you wanted to run a model without random intercepts by UnitID, and just random slopes, you'd need to run (0 + SamplingMonth | UnitID).

You mention that different FuelTypes would be expected to respond to seasonal changes differently. Right now you're modelling each individual stand as responding to seasonal changes differently. If you just wanted to model by FuelType, I think I would run something like this:

lmer(log(FMC) ~ Treatment + FuelType + (0 + SamplingMonth | FuelType) + (1 | UnitID), data=moisture)

Where you're specifying that you want a fixed effect of Treatment and FuelType, you want the effect of SamplingMonth to have different slopes in each FuelType, and you're modelling random intercepts for each UnitID. We don't model random intercepts for FuelType (thus the zero) because I assume from youre description that these are fixed categories and shouldn't be modeled as random representatives of an unknown distribution of Fuel Types. Hopefully that helps.

Update The model I posted above was a bit off. I ran the model below, which I believe is correct, on the actual data to no warnings - assuming MOISTURECONTENT is your outcome variable. Here we treat SamplingMonth as a random factor that has both a random intercept and random slopes of FuelType, which models the expected variance in seasonal changes across different FuelTypes.

summary(lmer(log(MOISTURECONTENT) ~ Treatment + FUELTYPE + (1 + FUELTYPE| SamplingMonth) + (1 |UNITID), data=moisture))

The results of this model indicate that there is substantially more variance in the moisture of the Litter Fueltype across months than the other FuelTypes. Note that this does treat months as being random factors, which is possibly not appropriate for your data if you expect there to be fixed effects of seasonal changes. But if you can assume these months are drawn from a random distribution of possible months within the season you tested, it should be sound. If you're concerned that that's not the case, you could potentially add a season factor (e.g. Autumn, Summer, Winter) and treat months within seasons as random factors.

If you run the code below, you'll be able to see the variation in the monthly effects by FuelType:

forestmodel <- lmer(log(MOISTURECONTENT) ~ Treatment + FUELTYPE + (1 + FUELTYPE| SamplingMonth) + (1 |UNITID), data=moisture)
require(lattice)
dotplot(ranef(forestmodel))

If it's of interest, adding an interaction term between the Treatment and FUELTYPE variables reveals that the Treatment seems to have most of its effect at the 1000hour and Litter Fueltypes. Not sure if that's relevant, but I found it interesting.

summary(lmer(log(MOISTURECONTENT) ~ Treatment*FUELTYPE + (1  + FUELTYPE| SamplingMonth) + (1 |UNITID), data=moisture))
$\endgroup$
  • $\begingroup$ @ Sean Thank you for your response and model suggestion. I ran the model, but it produced the error: ' Warning messages: 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00346073 (tol = 0.002, component 19) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?' $\endgroup$ – radex7 Oct 12 '14 at 20:06
  • $\begingroup$ <br/> The model executes fine when fixed effect FuelType is removed, or when I substitute FuelType with SamplingMonth.<br/>Thanks! $\endgroup$ – radex7 Oct 12 '14 at 20:20
  • $\begingroup$ It could be that there are too many parameters in your data for the relatively small number of groups (13). For instance, you might have FuelTypes with 1 or 2 forest stands in them. It might be easier if you could find a way to post/send the data to examine it more closely. $\endgroup$ – Sean Murphy Oct 13 '14 at 4:11
  • $\begingroup$ Sean, thank you so much! I am little embarrassed to admit, but per your suggestion I was formatting data to post online, and during the process I discovered a missing fuels sample (one row of data)! Upon re-entry I ran the model and I succeeded!!! Very happy camper right here! Just for potential future questions I posted a link to view and retrieve the data. Once again. Thank you so much! $\endgroup$ – radex7 Oct 13 '14 at 16:49
  • $\begingroup$ Thank for the nicely formatted link to your data! I ran the model I posted previously and got warnings, which made me realise that I'd mis-formatted part of it. To examine whether seasonal changes affect different fueltypes differently, you'd want to run (1 + FUELTYPE|SamplingMonth) rather than (0 + SamplingMonth | FuelType). The new one treats the months as the random factor where the fixed effects of FuelType might differ, which is what you want. The old one was treating Fueltypes as random factors where the effect of month might differ, and the model was having issues. I've updated my post. $\endgroup$ – Sean Murphy Oct 14 '14 at 1:21

Your Answer

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

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