Here's what I would do:
First, I would have a look here on how to specify the random term in your model1. I am not quite sure what you are trying to fit. There is also a lot of info on linear mixed effects models here on CV. Click on the lme4-nlme tag, which you also provided. It would also help if you could provide an example dataset, or at least the structure of your data.
Then, you most likely only need one model, which is presumably in the form of:
my_model <- lmer(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)
I specified the random effect to be + (1|site), because you said:
Study sites are included as the random effect in the model.
To get the ANOVA table you can either do:
library(car)
Anova(my_model)
or:
library(afex)
mixed(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)
or instead of running lmer() through the lme4 package, load the lmerTest package and run:
my_model <- lmer(carbon ~ species + landuse + species : landuse + (1|site), data = mydata)
anova(my_model)
This will give you the ANOVA table you probably need eventually. Make sure to have a look at those functions and their arguments (?Anova, ?mixed, ?lmerTest::anova).
I don't quite understand why would want to exclude species if the interaction is significant and run separate models for all species?!
However, if your main effects are not significant you could consider tossing them out and re-running the model with the interaction only. However, if one or both main effects are significant, I would keep them both in the model and report this together with a potential significant interaction.
In any case, if you have a significant interaction you should focus on interpreting the interaction and not the main effects since their interpretation could now be misleading. The interpretation of the interaction should start by visualizing it. You could do this for example using the emmip() function in the emmeans package:
library(emmeans)
emmip(my_model, landuse ~ species)
Regarding the adjustment of p-values, you only need to do that if you are following up with post-hoc tests.
This could be done with the emmeans() function (also from the emmeans package):
emmeans(my_model, pairwise ~ species : landuse)
Bendix Carstensen has a nice presentation in which he shows how this can be accomplished by splitting individual follow up time into small pieces and use a Poisson model for the rates thereby using the time scale (in your case age) as a covariate rather than part of the outcome. If you use R, the Epi package contains the tools you need.
Another option is to extend the Cox model to include time-dependent coefficients. See this vignette for the survival package in R. In that vignette it is also explained why it is not appropriate to include a covariate * age follow-up time interaction directly in the Cox model.
EDIT: Removed the bit about baseline age. That would of course not make sense with age as the time scale. Typed a bit too fast...
EDIT 2: Yes, the age * covariate interaction is an example of a time dependent coefficient and one of the ways for assessing the proportional hazards assumption.
As for the paper, it is possible that the critisism could apply. I can't find any info that the authors allowed for time-dependent coefficients. Including the age * covariate interaction is certainly possible, and R prints a warning message but estimates the coefficients. See the example below which is very inspired by the vignette mentioned above.
require(survival)
data(veteran)
## Use age as the time scale
## age * karno in the model
fit1 <- coxph(Surv(age, age + time, status == 1) ~ trt + celltype + karno + age:karno, data = veteran)
Warning message:
In coxph(Surv(age, age + time, status == 1) ~ trt + celltype + age * :
a variable appears on both the left and right sides of the formula
## karno as a time-dependent coefficient, simple interaction similar to fit1
fit2 <- coxph(Surv(age, age + time, status == 1) ~ trt + celltype + karno + tt(karno),
tt = function(x, t, ...) x * t, data = veteran)
## Coefficients from fit1
> fit1
Call:
coxph(formula = Surv(age, age + time, status == 1) ~ trt + celltype +
karno + age:karno, data = veteran)
coef exp(coef) se(coef) z p
trt 2.59e-01 1.30e+00 2.05e-01 1.26 0.2072
celltypesmallcell 8.32e-01 2.30e+00 2.73e-01 3.05 0.0023
celltypeadeno 1.17e+00 3.21e+00 2.95e-01 3.95 7.8e-05
celltypelarge 3.95e-01 1.48e+00 2.83e-01 1.39 0.1633
karno -3.33e-02 9.67e-01 1.07e-02 -3.10 0.0020
karno:age 4.75e-05 1.00e+00 1.77e-04 0.27 0.7890
Likelihood ratio test=62 on 6 df, p=1.75e-11
n= 137, number of events= 128
## Coefficients from fit2
> fit2
Call:
coxph(formula = Surv(age, age + time, status == 1) ~ trt + celltype +
karno + tt(karno), data = veteran, tt = function(x, t, ...) x *
t)
coef exp(coef) se(coef) z p
trt 1.70e-01 1.19e+00 2.04e-01 0.84 0.40333
celltypesmallcell 9.08e-01 2.48e+00 2.75e-01 3.30 0.00097
celltypeadeno 1.22e+00 3.40e+00 3.00e-01 4.08 4.5e-05
celltypelarge 4.04e-01 1.50e+00 2.83e-01 1.42 0.15422
karno -4.69e-02 9.54e-01 8.33e-03 -5.63 1.8e-08
tt(karno) 1.15e-04 1.00e+00 4.71e-05 2.45 0.01443
Likelihood ratio test=68.2 on 6 df, p=9.45e-13
n= 137, number of events= 128
Note the difference in the estimated coefficients. Reading the vignette again I see that section 4.2 deals with follow-up time * covariate interaction and not necessarily age * covariate interaction. But I think the example illustrates that results can be misleading when not properly allowing for time-dependent coefficients.
Hope this helps!
Best Answer
This is similar to having time-varying covariates in proportional-hazards survival modeling. The inherent assumption is that the association of the independent variable with outcome at a given time only depends on the current value of that variable at that time. The past history of that variable doesn't matter.
The interaction with time just allows the association of that independent variable (whatever its current value) with outcome to change over time since the start of the study. That isn't fundamentally different from any interaction term in a regression, as @Alexis notes in comments.
Whether that type of model makes sense in terms of the subject matter is another question. You might want to discuss that with those gave you this task.