# R gls() vs. SAS proc mixed with interaction: Why does R complain about a singular matrix when SAS does not?

I like to keep analyses all in SAS or all in R when I can help it and lately have been using R more and more, but there's one analysis that I do somewhat routinely that has given me trouble in R.

I have repeated measures data where I would like to fit the following model: $$Delta = Day + Group + Day\times Group$$ where $Delta$ is the change from baseline, $Day$ is the number of days from the beginning of the study, and $Group$ is the experimental group. I fit a variance-covariance matrix to account for repeated measures (for this example I'm using compound symmetry, but the difference is the same using other variance-covariance matrices). I have the data at the end of the post.

If I don't include the interaction, I can get the analysis to run as I want it to in both SAS and R. In SAS:

proc mixed data=df;
class group day id;
model delta = day group;
repeated day / subject=id type=cs;
lsmeans group / diff=all;
run;


In R:

library(nlme)
library(lsmeans)
fit.cs <- gls(Delta~Day+Group,
data=df,
corr=corCompSymm(,form=~1|ID))
anova(fit.cs,type="marginal")
lsmeans(fit.cs,pairwise~Group)


Obviously the results differ in terms of denominator DF, but I don't intend to start that discussion (unless that difference is causing the problem). When I add interaction in SAS, everything is fine:

proc mixed data=df;
class group day id;
model delta = day | group;
repeated day / subject=id type=cs;
run;


But when I do the same from R...

fit.cs <- gls(Delta~Day*Group,
data=df,
corr=corCompSymm(,form=~1|ID))
# Error in glsEstimate(object, control = control) :
#   computed "gls" fit is singular, rank 19


Why does R complain about the fit being singular but SAS doesn't?

Here are some fake data that are representative of data that I work with (from R):

df <- structure(list(Delta = c(-1.27, -0.34, 1.92, 0.45, 1.21, 0.43, -0.41, 0.16, -0.35,
1.49, -0.85, -0.86, 1.04, 0.49, 2.32, 0.13, -0.32, 0.5, 0.48, 1.21, -0.82, 0.93,
-0.58, 2.3, -0.9, 0.21, -0.72, 0.11, -0.28, -0.33, -0.7, -1.16, -0.23, -0.88, 0.97,
0.25, 0.8, 0.16, 0.63, -0.49, -0.63, -0.9, 1.1, -1.45, 0.38, -0.93, 0.4, 0.45, 0.48,
0.14, 1.02, -0.01, -1.98, 2.19, -1.53, -0.49, -1.57, -1.02, 1.09, 1.74, 0.54, -1.57,
-1.5, -0.48, 0.26, 0.2, -0.36, -1.05, -1.73, -0.77, -0.65, -1.07, -0.45, -0.14,
-0.56, 0.84, -2.66, -0.52, 1.44, 0.45, 0.24, -0.92), Day = structure(c(1L, 2L, 3L,
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L,
1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L,
3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 1L, 2L, 3L, 6L,
7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L, 1L, 2L, 3L, 6L, 7L), .Label = c("4",
"7", "10", "12", "14", "16", "28"), class = "factor"), Group = structure(c(1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("1", "2",
"3", "4"), class = "factor"), ID = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L,
4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 8L,
8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 11L,
12L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 15L, 15L, 15L,
15L, 15L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L),
.Label = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14",
"15", "16", "17", "18"), class = "factor")), .Names = c("Delta", "Day", "Group", "ID"),
class = "data.frame", row.names = c(NA, -82L))


Here's the same data for SAS:

DATA  df;
INPUT Delta Day Group $ID$ @@;
CARDS;
-1.27 4 1 1 -0.34 7 1 1 1.92 10 1 1 0.45 4 1 2 1.21 7 1 2 0.43 10 1 2
-0.41 4 1 3 0.16 7 1 3 -0.35 10 1 3 1.49 4 2 4 -0.85 7 2 4 -0.86 10 2
4 1.04 16 2 4 0.49 28 2 4 2.32 4 2 5 0.13 7 2 5 -0.32 10 2 5 0.5 16 2 5
0.48 28 2 5 1.21 4 2 6 -0.82 7 2 6 0.93 10 2 6 -0.58 16 2 6 2.3 28 2 6
-0.9 4 2 7 0.21 7 2 7 -0.72 10 2 7 0.11 16 2 7 -0.28 28 2 7 -0.33 4 2 8
-0.7 7 2 8 -1.16 10 2 8 -0.23 16 2 8 -0.88 4 3 9 0.97 7 3 9 0.25 10 3 9
0.8 16 3 9 0.16 28 3 9 0.63 4 3 10 -0.49 7 3 10 -0.63 10 3 10 -0.9 16 3 10
1.1 28 3 10 -1.45 4 3 11 0.38 7 3 11 -0.93 10 3 11 0.4 16 3 11 0.45 28 3 11
0.48 4 3 12 0.14 7 3 12 1.02 10 3 12 -0.01 16 3 12 -1.98 28 3 12 2.19 4 3 13
-1.53 7 3 13 -0.49 10 3 13 -1.57 16 3 13 -1.02 28 3 13 1.09 4 4 14 1.74 7 4 14
0.54 10 4 14 -1.57 16 4 14 -1.5 4 4 15 -0.48 7 4 15 0.26 10 4 15 0.2 16 4 15
-0.36 28 4 15 -1.05 4 4 16 -1.73 7 4 16 -0.77 10 4 16 -0.65 16 4 16 -1.07 28 4 16
-0.45 4 4 17 -0.14 7 4 17 -0.56 10 4 17 0.84 16 4 17 -2.66 28 4 17 -0.52 4 4 18
1.44 7 4 18 0.45 10 4 18 0.24 16 4 18 -0.92 28 4 18
;
RUN;

• I ran into the same problem when performing an exploratory factor analysis a year ago! SAS would run like everything was fine (even though some variables were exact linear combinations of other variables included), but R would (appropriately) throw an error Mar 26, 2014 at 15:11
• In the R code, you have delta, and Deltaa and then grp and Group. I can't correct this because edits must be 6 characters. Mar 26, 2014 at 16:15

## 2 Answers

This is only a partial answer so far. It seems the main culprit is the fact that not all of the groups have data at all of the timepoints, which is what causes the matrix to be singular when including the interaction term, so I can get similar results by running:

fit.cs <- gls(Delta~Day*Group,
data=df[df\$Day %in% c(4,7,10),],
corr=corCompSymm(,form=~1|ID))


and

proc mixed data=df;
where Day in (4,7,10);
class group day id;
model delta = day | group;
repeated day / subject=id type=cs;
run;


The question about what SAS does to get around the lack of information at particular time points still remains, but having discovered the main trigger will help me move toward a happy resolution.

I got a similar error message and finally discovered that I had too much missing data in the y-variable of my dataset for certain factor categories. When I replaced many of my NA's with numbers (luckily in my case those NA's were supposed to be zeros), the model ran fine.