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;