As others have pointed out, the problem is singularity that the estimated variance of random intercepts is at its boundary zero. In such cases, lme4 prints messages but nlme does not. However, you can use nlme::interval() to assess boundary issues. In fact, your model will produce the following error when finding confidence intervals.
library(nlme)
library(lme4)
library(ggplot2)
library(dplyr)
Data <- read.table(header = TRUE, stringsAsFactors = TRUE, text = "
UID MID Seq Score Condition
1 1 1 1 90 B
2 1 1 2 85 A
3 1 1 3 75 B
4 2 4 1 75 A
5 2 4 2 95 B
6 2 4 3 85 A
7 3 6 1 60 A
8 3 6 2 82 B
9 3 6 3 85 A
10 4 3 1 60 A
11 4 3 2 70 B
12 4 3 3 75 A
13 5 2 1 85 B
14 5 2 2 85 A
15 5 2 3 85 B
16 6 5 1 90 B
17 6 5 2 95 A
18 6 5 3 100 B
19 7 2 1 90 B
20 7 2 2 70 A
21 7 2 3 50 B
22 8 1 1 70 B
23 8 1 2 75 A
24 8 1 3 80 B
25 9 3 1 90 A
26 9 3 2 30 B
27 9 3 3 90 A
28 10 6 1 50 A
29 10 6 2 85 B
30 10 6 3 92 A
31 11 4 1 50 A
32 11 4 2 85 B
33 11 4 3 92 A
34 12 5 1 65 B
35 12 5 2 50 A
36 12 5 3 90 B
37 13 4 1 65 A
38 13 4 2 70 B
39 13 4 3 80 A
40 14 2 1 60 B
41 14 2 2 100 A
42 14 2 3 80 B
")
summary(Model <- lme(Score ~ Condition, random = ~ 1 | UID, data = Data))
intervals(Model)
"Error in intervals.lme(Model) :
cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance
Consider which = 'fixed'"
It means that the estimated variance is not positive, which is in line with a zero-variance case. The reason for singularity in this case is not incorrect syntax but the data characteristics. Before modeling, we need to explore the data in tables and plots.
sd(Data$Score)
"15.79279"
Data %>%
group_by(UID) %>%
summarize(Umean = mean(Score), Usd = sd(Score)) %>%
summarize(REsd = sd(Umean), Usdm = mean(Usd))
" UID Umean Usd
<int> <dbl> <dbl>
1 1 83.3 7.64
2 2 85 10
3 3 75.7 13.7
4 4 68.3 7.64
5 5 85 0
6 6 95 5
7 7 70 20
8 8 75 5
9 9 70 34.6
10 10 75.7 22.5
11 11 75.7 22.5
12 12 68.3 20.2
13 13 71.7 7.64
14 14 80 20
REsd 7.78
Usdm 14.0"
Data %>%
group_by(UID) %>%
mutate(DM = Score - mean(Score)) %>%
ungroup() %>%
summarize(DMsd = sd(DM))
"13.8"
The standard deviation of Score is 15.8. After removing means by user, the standard deviation is on average 14.0, or more accurately 13.8, which does not reduce much. The standard deviation of user-specific means is 7.78, much smaller than the standard deviation of demeaned scores and will shrink further in random-effects estimators. Recall that random effects are user-specific sequence-average deviation from the ground mean that should make the within-subject variation smaller than the total variation. If within-subject variation equals the total variation, there are no random effects.
Data <- Data %>%
group_by(UID) %>%
mutate(Group = if_else(
Condition[1] == "B", "BAB", "ABA") %>%
factor(levels = c("ABA", "BAB"))) %>%
ungroup()
ggplot(
data = Data,
aes(x = factor(Seq), y = Score, group = UID, color = Group) ) +
geom_smooth(aes(
ymin = ifelse(after_stat(ymin) < 0, 0, after_stat(ymin)),
ymax = ifelse(after_stat(ymax) > 100, 100, after_stat(ymax)),
group = Group, fill = Group),
method = "loess", span = 1, linetype = 0, show.legend = F) +
stat_smooth(
geom = "line", aes(group = Group),
method = "loess", span = 1, linewidth = 1, alpha = 0.7) +
geom_line(linewidth = 0.2, alpha = 0.4, show.legend = F) +
scale_fill_manual(values = c("#2197DF", "#D14343")) +
scale_color_manual(values = c("#2197DF", "#D14343")) +
scale_x_discrete(expand = c(0, 0.1, 0, 0.1)) +
labs(x = "Sequence") +
theme_minimal() +
theme(
legend.title = element_blank(),
legend.position = c(0.6, 0.2),
panel.grid = element_line(colour = "grey90", linewidth = 0.2)
)
We can also view such data characteristics from a spaghetti plot. It shows no patterns that user-specific mean values deviate from the ground mean. In contrast, the variation within each user is so great that it is comparable with total variation. However, the comparison should not be made based on Condition but the first Condition each user experienced that sufficiently defines the two distinct experimental sequences, which is coded into Group. Groups ABA and BAB appear to have different means at three sequential stages that show linear effects on Score. Starting with B, the delay is longer than that starting with A but remains constant in additional stages, unlike those starting with A whose delay increases in the second and third stages. So we can use Seq as numeric instead of a factor and should have interaction between Seq and Group to allow different intercepts and slopes of Seq by Group.
summary(Model <- lmer(Score ~ 1 + (1 | UID), data = Data))
summary(Model <- lmer(Score ~ Condition + (1 | UID), data = Data))
summary(Model <- lmer(Score ~ factor(Seq) + (1 | UID), data = Data))
summary(Model <- lmer(Score ~ Group + (1 | UID), data = Data))
summary(Model <- lmer(Score ~ Group * Seq + (1 | UID), data = Data))
"boundary (singular) fit: see help('isSingular')
REML criterion at convergence: 324.2
Groups Name Variance Std.Dev.
UID (Intercept) 0.0 0.00
Residual 220.4 14.85
Number of obs: 42, groups: UID, 14
Estimate Std. Error t value
(Intercept) 53.286 8.571 6.217
GroupBAB 24.810 12.122 2.047
Seq 10.643 3.968 2.682
GroupBAB:Seq -9.929 5.611 -1.769"
Because the within-subject variation is as large as the total variation, any attempt to estimate random intercepts is challenging. lme4::lmer() shows boundary issues and zero variance in all specifications with random intercepts by UID.
summary(Model <- gls(
Score ~ Group * Seq, correlation = corCompSymm(form = ~ 1 | UID), data = Data))
intervals(Model, which = "var-cov")
" Correlation structure:
lower est. upper
Rho -0.3180153 -0.1002069 0.2332579"
summary(Model <- gls(
Score ~ Group * Seq, correlation = corAR1(form = ~ Seq | UID), data = Data))
intervals(Model, which = "var-cov")
" Correlation structure:
lower est. upper
Phi -0.3811611 -0.05991643 0.2742381"
Without using random effects, we can still incorporate clustered errors by specifying error correlation patterns in generalized least squares. The above summary shows that the correlation is no different from zero in either a compound symmetry structure or an autoregressive process. Note that neither estimates nor confidence intervals of correlation coefficients are at boundary -1 or 1, which implies stable results.
summary(Model <- gls(
Score ~ Group * Seq, data = Data))
" AIC BIC logLik
334.2345 342.4224 -162.1173
Value Std.Error t-value p-value
(Intercept) 53.28571 8.571399 6.216688 0.0000
GroupBAB 24.80952 12.121789 2.046688 0.0477
Seq 10.64286 3.967787 2.682316 0.0108
GroupBAB:Seq -9.92857 5.611298 -1.769390 0.0849"
summary(Model <- gls(
Score ~ Group * Seq, weights = varIdent(form = ~ 1 | Seq), data = Data))
"Generalized least squares fit by REML
AIC BIC logLik
335.685 347.1481 -160.8425
Variance function:
Structure: Different standard deviations per stratum
Formula: ~1 | Seq
Parameter estimates:
1 2 3
1.0000000 1.3269816 0.8687116
Value Std.Error t-value p-value
(Intercept) 53.36942 7.781308 6.858670 0.0000
GroupBAB 24.67001 11.004431 2.241825 0.0309
Seq 10.67273 3.395007 3.143656 0.0032
GroupBAB:Seq -9.97837 4.801265 -2.078279 0.0445
Residual standard error: 13.5885
Degrees of freedom: 42 total; 38 residual"
The spaghetti plot shows that both ABA and BAB groups have similar spread in each stage, and both groups have slightly larger spread in the second stage. Therefore, we can use weighted least squares to allow error variance to vary by stage. As shown above, compared to OLS and random effects, WLS results in smaller standard errors and p values. However, both AIC and BIC prefer OLS.
In summary, not all longitudinal data require random-effects or fixed-effect estimators. When the between-subject deviation is small, within-subject variation is comparable with total variation, random effects approach zero variance, and simpler models may give better results.
Seqas a covariate I can have an idea about it. However, with or without that, the procedure fails to fit the random effects. $\endgroup$