2
$\begingroup$

I am using the lme4 package to create a linear mixed effects model that accounts for repeated measures.

Data

colony_area_and_cover_code_data <- structure(list(TimeStep = structure(c(3L, 3L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 
7L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 
9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 
10L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 
12L, 12L, 12L, 12L, 12L, 12L, 12L), .Label = c("1", "2", "3", 
"4", "5", "6", "7", "8", "9", "10", "11", "12"), class = "factor"), 
    Site_long = c("Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Waikiki", "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Waikiki", "Waikiki"
    ), Shelter = c("High", "Low", "High", "High", "Low", "Low", 
    "High", "High", "Low", "Low", "High", "High", "Low", "Low", 
    "High", "High", "High", "Low", "Low", "High", "High", "Low", 
    "Low", "Low", "High", "High", "High", "Low", "Low", "High", 
    "High", "Low", "Low", "Low", "High", "High", "High", "Low", 
    "Low", "High", "High", "Low", "Low", "Low", "High", "High", 
    "High", "Low", "Low", "High", "High", "Low", "Low", "Low", 
    "High", "High", "High", "Low", "Low", "High", "High", "Low", 
    "Low", "Low", "High", "High", "High", "Low", "Low", "High", 
    "High", "Low", "Low", "Low", "High", "High", "High", "Low", 
    "Low", "High", "High", "Low", "Low", "Low", "High", "High", 
    "Low", "Low"), `Module #` = structure(c(10L, 9L, 8L, 10L, 
    9L, 11L, 3L, 5L, 4L, 6L, 8L, 10L, 9L, 11L, 1L, 3L, 5L, 4L, 
    6L, 8L, 10L, 7L, 9L, 11L, 1L, 3L, 5L, 2L, 6L, 8L, 10L, 7L, 
    9L, 11L, 1L, 3L, 5L, 2L, 6L, 8L, 10L, 7L, 9L, 11L, 1L, 3L, 
    5L, 2L, 6L, 8L, 10L, 7L, 9L, 11L, 1L, 3L, 5L, 2L, 6L, 8L, 
    10L, 7L, 9L, 11L, 1L, 3L, 5L, 2L, 6L, 8L, 10L, 7L, 9L, 11L, 
    1L, 3L, 5L, 2L, 6L, 8L, 10L, 7L, 9L, 11L, 1L, 3L, 2L, 6L), .Label = c("111", 
    "112", "113", "114", "115", "116", "212", "213", "214", "215", 
    "216", "211"), class = "factor"), total_area = c(10.4730161798765, 
    0.470417310950379, 1.12170962265947, 15.1461766957539, 0.84071836278275, 
    1.3177670552551, 0.653559512764809, 0.370705324910543, 1.14688997513326, 
    0.581990145895235, 2.15410733439833, 14.3418757404116, 1.74013649762257, 
    4.61887925562605, 0.0215112373734017, 1.12357255063702, 0.172065091620736, 
    1.88495559215388, 0.84071836278275, 3.92523172214192, 18.822153496595, 
    0.350263995996734, 0.320718840997778, 6.3469169931263, 0.0904658706496261, 
    1.50050498439191, 2.47944359056696, 0.0500581295118801, 1.43055270954443, 
    7.03223926207256, 28.3286133694053, 0.210542714934101, 0.464586529499474, 
    5.16045375059151, 0.788854045723908, 1.36454209678304, 1.08502177635188, 
    1.36439318798517, 1.24222016364484, 13.6424060473263, 34.5765397419597, 
    1.87898302475959, 2.09734513325649, 3.28787654577659, 2.45098999462753, 
    0.180931741299252, 0.46991974190878, 2.13621743874264, 1.33988663867242, 
    16.1393724864509, 28.4216296967568, 3.14191620673555, 3.02078382765752, 
    0.569662911350172, 2.6737313910935, 0.461171195560169, 0.282854187854266, 
    3.04181900058373, 0.399158921101302, 13.6005080948471, 16.4004774522405, 
    3.89003406709318, 0.44919611587226, 0.895698013383031, 3.47412025908296, 
    0.370705324910543, 0.233200331381288, 2.8707558074787, 1.40500001668743, 
    10.6809732323847, 9.62324008176214, 14.1008970417224, 4.82127391024982, 
    3.83762847457455, 0.746640117023492, 0.0430224747468034, 
    0.0904658706496261, 2.88295085172956, 0.944133971239514, 
    8.82540502985815, 9.94146284472563, 9.86581183076813, 7.70728361307783, 
    0.891074937564863, 0.653559512764809, 0.192792590243542, 
    1.42572751536521, 0.346431271895304), Season = c("winter", 
    "winter", "spring", "spring", "spring", "spring", "spring", 
    "spring", "spring", "spring", "summer", "summer", "summer", 
    "summer", "summer", "summer", "summer", "summer", "summer", 
    "fall", "fall", "fall", "fall", "fall", "fall", "fall", "fall", 
    "fall", "fall", "winter", "winter", "winter", "winter", "winter", 
    "winter", "winter", "winter", "winter", "winter", "spring", 
    "spring", "spring", "spring", "spring", "spring", "spring", 
    "spring", "spring", "spring", "summer", "summer", "summer", 
    "summer", "summer", "summer", "summer", "summer", "summer", 
    "summer", "fall", "fall", "fall", "fall", "fall", "fall", 
    "fall", "fall", "fall", "fall", "winter", "winter", "winter", 
    "winter", "winter", "winter", "winter", "winter", "winter", 
    "winter", "spring", "spring", "spring", "spring", "spring", 
    "spring", "spring", "spring", "spring"), Year = c(17, 17, 
    18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 
    18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 18, 
    18, 18, 18, 18, 18, 18, 18, 19, 19, 19, 19, 19, 19, 19, 19, 
    19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 
    19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 19, 
    19, 19, 20, 20, 20, 20, 20, 20, 20, 20, 20), mean_cover_code = c(1.03703703703704, 
    1.75, 1.125, 1.11428571428571, 1.09090909090909, 1.13793103448276, 
    1.2962962962963, 1.2258064516129, 1.5, 1.29411764705882, 
    1.2, 1.05333333333333, 1, 1.225, 1.4, 1.65217391304348, 1.75675675675676, 
    1.71428571428571, 1.58333333333333, 1, 1.00934579439252, 
    1.64705882352941, 1.75, 1.22857142857143, 1.64, 1.66666666666667, 
    1.2, 1.04545454545455, 1.26785714285714, 1.38235294117647, 
    1.06766917293233, 1.85714285714286, 1.57142857142857, 2, 
    1.37142857142857, 1.25806451612903, 1.14285714285714, 1.33333333333333, 
    1.28301886792453, 1.425, 1.44067796610169, 1.20833333333333, 
    1.2962962962963, 1.1219512195122, 1.425, 1.29268292682927, 
    1.45652173913043, 1.22058823529412, 1.26086956521739, 1.27659574468085, 
    1.14285714285714, 1.10526315789474, 1.63636363636364, 1.83870967741935, 
    1.40816326530612, 1.92682926829268, 1.18181818181818, 1.7962962962963, 
    1.35714285714286, 2, 1.14492753623188, 1.70588235294118, 
    1.95, 1.13333333333333, 1.45945945945946, 1.85714285714286, 
    1.36585365853659, 1.42222222222222, 1.44897959183673, 1.27906976744186, 
    1.61538461538462, 2, 1.56410256410256, 1.59375, 1.88, 2.24242424242424, 
    2.18181818181818, 2.25, 2.08333333333333, 1.1875, 1.12903225806452, 
    1.125, 1.32142857142857, 1.0625, 1.45454545454545, 1.68421052631579, 
    1.4, 1.16666666666667), Date_new = structure(c(17484, 17484, 
    17575, 17575, 17575, 17575, 17600.8, 17594, 17601, 17586, 
    17687, 17680, 17687, 17683.5, 17684, 17675, 17677, 17682, 
    17675, 17783, 17783, 17783, 17783, 17783, 17780, 17782, 17773, 
    17773, 17782, 17869, 17862, 17862, 17862, 17862, 17860, 17859, 
    17867, 17867, 17860, 17960, 17967, 17953, 17960, 17953, 17964, 
    17965, 17957, 17957, 17964, 18044, 18037, 18037, 18044, 18044, 
    18046, 18039, 18048, 18047, 18046, 18142, 18142, 18135, 18135, 
    18135, 18144, 18139, 18140, 18140, 18144, 18233, 18233, 18226, 
    18226, 18226, 18231, 18223, 18230, 18230, 18231, 18333, 18333, 
    18333, 18331, 18331, 18315, 18321, 18336, 18315), class = "Date")), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -88L), groups = structure(list(
    TimeStep = structure(c(3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
    5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 
    9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 12L, 12L, 
    12L, 12L), .Label = c("1", "2", "3", "4", "5", "6", "7", 
    "8", "9", "10", "11", "12"), class = "factor"), Site_long = c("Hanauma Bay", 
    "Hanauma Bay", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Hanauma Bay", "Hanauma Bay", 
    "Waikiki", "Waikiki", "Hanauma Bay", "Hanauma Bay", "Waikiki", 
    "Waikiki", "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", 
    "Hanauma Bay", "Hanauma Bay", "Waikiki", "Waikiki", "Hanauma Bay", 
    "Hanauma Bay", "Waikiki", "Waikiki", "Hanauma Bay", "Hanauma Bay", 
    "Waikiki", "Waikiki"), Shelter = c("High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low", "High", "Low", "High", "Low", "High", 
    "Low", "High", "Low"), .rows = structure(list(1L, 2L, 3:4, 
        5:6, 7:8, 9:10, 11:12, 13:14, 15:17, 18:19, 20:21, 22:24, 
        25:27, 28:29, 30:31, 32:34, 35:37, 38:39, 40:41, 42:44, 
        45:47, 48:49, 50:51, 52:54, 55:57, 58:59, 60:61, 62:64, 
        65:67, 68:69, 70:71, 72:74, 75:77, 78:79, 80:81, 82:84, 
        85:86, 87:88), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -38L), .drop = TRUE))

In this analysis, I am looking to see if two variables are correlated, namely if algal cover (mean_cover_code) is a significant predictor of coral cover (total_area) within 11 experimental modules that were sampled repeatedly over time where 6 are at one site (Waikiki) and 5 are at another (Hanauma Bay). In order to account for repeated measures, I created random effects for these experimental modules (1|module_colony_area_vs_cover_code) and time with nested season inside year (1|Season/Year).

First model with mistakenly switched predictor and response

# Variable vector of module number to account for repeated measures
module_colony_area_vs_cover_code <- colony_area_and_cover_code_data$`Module #`

# Mistake model
col_area_vs_cover_code_lmer_only <- lmer(mean_cover_code ~ total_area + (1|module_colony_area_vs_cover_code) + (1|Year/Season), data = colony_area_and_cover_code_data, na.action = "na.fail")
summary(col_area_vs_cover_code_lmer_only)

Originally when I did this analysis, I made the mistake of making coral cover a predictor of algal cover on accident. This model (as the output code is shown above) showed a significant relationship with coral cover as the predictor (total_area) and algal cover as the response (mean_cover_code).

Once I realized I made this mistake, I switched algal cover to be the predictor and coral cover to be the response. However I was shocked to find out when I did this, there was a gigantic change in the p-value generated by the lmer fit. Whereas the mistaken model had a p-value equal to $0.0202$ (see mistake model output), the corrected model with the predictor and the response reversed was $p = 0.887$. See the below script for the output:

Second model with correct predictor and response

# Variable vector of module number to account for repeated measures
module_colony_area_vs_cover_code <- colony_area_and_cover_code_data$`Module #`

# Correct model
col_area_vs_cover_code_lmer_only <- lmer(total_area ~ mean_cover_code + (1|module_colony_area_vs_cover_code) + (1|Year/Season), data = colony_area_and_cover_code_data, na.action = "na.fail")
summary(col_area_vs_cover_code_lmer_only)

I was very confused as to how the significance of a lmer linear mixed effects model could change so drastically with swapping the predictor and response variables. Especially considering that looking at the plot of the raw data, there does appear to be a correlation despite the model results.

Figure Code

col_area_vs_cover_code_plot <- ggplot(data = colony_area_and_cover_code_data, aes(x = mean_cover_code, y = total_area, fill = interaction(Site_long, Shelter), shape = interaction(Site_long, Shelter))) +
  geom_point(aes(size = 5)) +
  stat_smooth(method=stats::lm, aes(fill = NULL, shape = NULL)) +
  scale_shape_manual(name = 'Site x Shelter', values = c(21, 24, 21, 24), labels = c("Hanauma Bay - Low", "Waikiki - Low", "Hanauma Bay - High", "Waikiki - High")) +
  scale_fill_manual(name = "Site x Shelter", values = c("black", "black", "white", "white"), labels = c("Hanauma Bay - Low", "Waikiki - Low", "Hanauma Bay - High", "Waikiki - High")) +
  guides(size = FALSE, linetype = FALSE,  shape = guide_legend(override.aes = list(size = 4))) +
  labs(x = "Algal overgrowth (1-4)", y = expression(paste("Total coral cover (cm"^"2",")"))) +
  #coord_cartesian(xlim = c(1, 1.7)) +
  theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.title = element_text(size = 16), axis.text.y = element_text(angle = 90), axis.text = element_text(size = 16), legend.text = element_text(size = rel(1.5)), legend.title = element_text(size = rel(1.5)), legend.position = "none", plot.title = element_text(size = 20, hjust = 0.5, vjust = -1.5))

enter image description here enter image description here

I wanted to see if anyone could explain how this is happening or if something is wrong with the analysis and/or figure code that is making a correlation completely change when swapping the predictor and response variable. Especially because my understanding that when two variables are correlated with one another, that correlation should exist regardless of which is the predictor and which is the response.

$\endgroup$
5
  • $\begingroup$ GLMMs are an extensions of GLMs, which with a Gaussian distribution family and an identity link are the same as ordinary least-squares regression models. It's easier to understand why this happens if you step back and revisit OLS regression. What exactly is minimized in OLS regression and how could that be important to your question? $\endgroup$
    – Roland
    Commented Sep 15, 2023 at 5:38
  • 1
    $\begingroup$ @Roland: swapping $X$ and $Y$ in OLS regression doesn't change the $p$-value for association, so I don't see why it's going to be very illuminating about why the $p$-value does change for a mixed model $\endgroup$ Commented Sep 15, 2023 at 6:20
  • $\begingroup$ @ThomasLumley That's correct but the actual fit changes. In a GLMM, both the estimates of fixed- and random-effects parameters will change (OP's second model even gives a singular fit because a random effect is estimated as zero). Why would you expect to have identical p-values for the fixed-effect slope then? $\endgroup$
    – Roland
    Commented Sep 15, 2023 at 6:40
  • 1
    $\begingroup$ I wouldn't, but that's because this is not like OLS rather than because it is like OLS. $\endgroup$ Commented Sep 15, 2023 at 6:47
  • $\begingroup$ With the correct model, doesn't the coefficient for mean_cover_code change sign as well? In the second model the estimate for mean_cover_code is positive (I get 0.2154) which contradicts the picture and I guess also your intuition? I'm wondering whether this is more interesting than what some p-value. $\endgroup$
    – dipetkov
    Commented Sep 17, 2023 at 6:47

2 Answers 2

6
$\begingroup$

From the point of view of the fixed part of a linear mixed model, the random part of the model controls how different points are weighted in the regression. When you have the model the correct way around, you have a fairly large random intercept standard deviation for Module # (it's about a third larger than the residual standard devaition).

That implies quite strong correlation within module, so that pairs of observations in the same module will be given more weight in the regression if they differ in mean_cover_code and less weight if they have the same mean_cover_code than in a linear regression. On the other hand, the Year and Season standard deviations are pretty small, so Year and Season don't affect the weighting all that much.

When the model is fitted backwards, the standard deviation corresponding to Module # is basically zero. The only variance component of any real size is the Season in Year component, and its standard deviation is only about half the size of the residual. So there's a bit of upweighting/downweighting by Season, but basically none by Module.

Comparing to ordinary linear regression, then, you might expect the correct model (with big variance components) to disagree with an ordinary linear model as to the regression slope (and it does). In the backwards model (without much up/downweighting) you might expect the mixed model to agree pretty well with the ordinary linear regression slope (and it does).

Somehow, this difference in weighting is explaining the difference in fitted regression. Given that the association in your graph looks to be driven by the difference between the black circles and the other points, that's plausible. I can imagine that the variance part of the mixed model leads to downweighting the comparison between black circles and other points to some extent, and this make the association go away.

(There will also be differences in p-values because the standard error estimates are different, but the differences in the sign of the slope are the thing that really needs explaining)

$\endgroup$
1
  • $\begingroup$ Thank you very much for your detailed response. I wanted to provide an updated figure to you that labels all the data points with their corresponding Module # to potentially illustrate your point (see second figure in edited post). In your opinion, which modules are most attributing to this issue? I have been trying to look for leverage points in this analysis to understand how I could explain what is going on in a biological context. The hypothesis is that we expect less coral cover with greater algal overgrowth. $\endgroup$ Commented Sep 16, 2023 at 21:43
1
$\begingroup$

The dataset is a lot more interesting and rich that the "correct" model suggests. Here is a small-multiples plot which highlights the differences between sites, seasons and shelter types.

enter image description here

The two sites, Hanauma Bay and Waikiki, seem to operate under different regimes. I'm not sure it makes sense to analyze the two sites together.

So let's focus on Hanauma Bay, which exhibits more variability in total coral cover. Now I also show the module IDs; for readability abbreviate "212" as "2", "213" as "3", etc.

enter image description here

It's clear that two modules ("213" and "215" abbreviated as "3" and "5") provide high shelter and the rest — low shelter.

Intuitively, there is less variability in algal growth in the two high-shelter modules. In other words, the negative correlation between algal overgrowth and coral cover is driven to a large extent by the omission of shelter type. Perhaps you can make a scientific argument this is a reasonable thing to do?

Finally, you can consider Generalized Least Squares (GSL) as an alternative to the Linear Mixed Model (LMM). Here is how I specify a model for total area as a function of mean cover, with Module-specific variances and continuous autocorrelation on time.

library("nlme")
fit.gsl <- gls(
  total_area ~ mean_cover * I(Season == "spring"),
  weights = varIdent(form = ~ 1 | Module),
  correlation = corCAR1(form = ~ Date | Module),
  data = data2
)
summary(fit.gsl)

#> Coefficients:
#>                                           Value Std.Error    t-value p-value
#> (Intercept)                            5.387768  2.935519  1.8353717  0.0737
#> mean_cover                            -1.437376  1.862865 -0.7715939  0.4448
#> I(Season == "spring")TRUE            -26.769467 10.665838 -2.5098327  0.0161
#> mean_cover:I(Season == "spring")TRUE  22.429773  8.948023  2.5066736  0.0162

I've added an interaction with Season == "spring" since the spring pattern is qualitatively different than the other seasonal patterns. The p-value for algal overgrowth is 0.4 but the graphs already showed that the variance within a module is high: the residual std. error is 5.6 while the cover effect is only -1.4.

At least the coefficient estimate is in the "right" direction (negative) which is actually not true in your "correct" model where the effect of cover is 0.2154.

PS: The GLS results are about the same with a corAR1(form = ~ 1 | Module) correlation. This choice is also reasonable (more reasonable?) than the corCAR1 structure: Even though the earliest observation is taken at a different time for each module, once a module is "selected" for the study, there is exactly one observation for each (Season, Year) combination until Spring 2020.

In any case, the important difference between the LMM and the GLS models is how they model the within-module variance. The GSL model allows each module to have its own variance:

#> Variance function:
#>  Structure: Different standard deviations per stratum
#>  Formula: ~1 | Module 
#>  Parameter estimates:
#>       212       213       214       215       216 
#> 1.0000000 1.2197720 0.4324221 2.9470512 0.3749878 

These variance estimates agree well with the second plot above which shows that the two high-shelter modules (213 and 215) have higher variance in coral cover than the three low-shelter modules (212, 214, 216).

The LMM doesn't quite capture this situation because it assumes the random module effects have the same variance:

ranef(fit.lmm)$Module
#>     (Intercept)
#> 212  -2.1375317
#> 213   0.4338779
#> 214  -5.3925046
#> 215  10.1676394
#> 216  -3.0714810
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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