R gls() vs. SAS proc mixed con interacción: ¿Por qué R se queja de una matriz singular cuando SAS no lo hace?

Question

Me gusta mantener los análisis todo en SAS o todo en R cuando puedo evitarlo y últimamente he estado usando R cada vez más, pero hay un análisis que hago de forma algo rutinaria que me ha dado problemas en R.

Tengo datos de medidas repetidas en los que me gustaría ajustar el siguiente modelo: $$Delta = Day + Group + Day\times Group$$ donde $Delta$ es el cambio respecto a la línea de base, $Day$ es el número de días desde el inicio del estudio, y $Group$ es el grupo experimental. Ajusté una matriz de varianza-covarianza para tener en cuenta las medidas repetidas (para este ejemplo estoy utilizando la simetría compuesta, pero la diferencia es la misma utilizando otras matrices de varianza-covarianza). Tengo los datos al final del post.

Si no incluyo la interacción, puedo conseguir que el análisis se ejecute como quiero tanto en SAS como en R. En SAS:

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

En 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)

Evidentemente, los resultados difieren en cuanto a la DF del denominador, pero no pretendo iniciar esa discusión (a menos que esa diferencia sea la causa del problema). Cuando añado la interacción en SAS, todo va bien:

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

Pero cuando hago lo mismo desde 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

¿Por qué R se queja de que el ajuste es singular pero SAS no?

Aquí hay algunos datos falsos que son representativos de los datos con los que trabajo (de 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))

Aquí están los mismos datos para 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;

Answer 1

Esta es sólo una respuesta parcial hasta ahora. Parece que el principal culpable es el hecho de que no todos los grupos tienen datos en todos los puntos de tiempo, que es lo que hace que la matriz sea singular al incluir el término de interacción, por lo que puedo obtener resultados similares corriendo:

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

y

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;

La pregunta sobre qué hace SAS para sortear la falta de información en determinados momentos sigue en pie, pero haber descubierto el principal desencadenante me ayudará a avanzar hacia una feliz resolución.

Answer 2

Recibí un mensaje de error similar y finalmente descubrí que tenía demasiados datos perdidos en la variable y de mi conjunto de datos para ciertas categorías de factores. Cuando sustituí muchos de mis NA por números (por suerte, en mi caso esos NA debían ser ceros), el modelo funcionó bien.