Controlling for an effect by adding it as covariate in R

Question

I have a dataset like this:

individual1	individual2	Bray_Curtis	Dyad_type	proximity	matriline	id1	id2
Aapi	App	0.47	1	0.14	1	1	2
Aapi	Eis	0.60	2	0.03	0	1	4
Potj	Popp	0.50	1	0.11	1	3	5
Aapi	Potj	0.27	3	0.15	0	1	3
Ree	Eis	0.63	4	0.66	0	6	4

I want to look at a possibile effect of the dyad type covariate on diet similarity (bray curtis indices). I have 4 levels of 'dyad_type', defined depending on the relationship between individual1 and individual2. The variable 'proximity' is the percent of times the individuals in the dyad were in proximity from each other when the other one was eating, and 'matriline' means that the 2 individuals are part (1) or not (0) of the same matriline.

I know that dyad_type and matriline affect proximity:

lmm1 <- lmer(proximity~dyad_type + matriline + (1|id1) + (1|id2))

Analysis of Deviance Table (Type III Wald chisquare tests)

Response: proximity
             Chisq Df Pr(>Chisq)    
(Intercept) 58.733  1  1.806e-14 ***
dyad_type   16.750  3  0.0007955 ***
matriline   19.777  1  8.702e-06 ***

these results make sense, as we expect some individuals to be more tolerant to each other while feeding, hence more in close spatial proximity, and members of the same matriline are expected to be generally more in close proximity.

HOWEVER, the aim of my analysis is to investigate if different dyad types and if individuals from the same matriline have more or less similar diets (described by bray-curtis indices), and I want to control for the effect from proximity.

Here is my model that seems to explain my data better (based on model comparisons with more or less interaction terms):

lmm2 <- lmer(formula = bray_curtis ~ dyad_type + proximity + matriline + 
                        proximity*dyad_type + (1|id1) + (1|id2)

Am I doing the right thing here? I'm not sure if this is the right way to control for proximity in my results.

Thank you!

Answer 1

That looks quite good.
Just a few notes:

First, you should always visualize your data. And if you deal with interactions, consider using some of the interactive visualization tools, see e.g. here or here.

Second, make sure that the variables dyad_type, matriline, id1, and id2 are coded as factors.

Third, in the formula you posted for lmm2, if you use the interaction term proximity*dyad_type, you can drop the proximity and the dyad_type term, both formulae would result in the same model.

Fourth, the easiest way to control for proximity would be to only use the term proximity. But if your formula works better, i.e. the data fit is much better, the interaction term is justified. It simply means that the proximity slope differs considerably for different values of dyad_type. If that corresponds to your domain knowledge, you should leave it that way. Just keep in mind, that interactions make your model considerably more complex and interpretation more difficult. Again, checking your visualizations should help you make a decision.

Fifth, it is interesting that you use an interaction proximity*dyad_type but no interaction proximity*matriline. But if you have ruled it out because it didn't significantly improve the data fit, great.