4 votos

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

Me gusta mantener todos los análisis en SAS o todos en R cuando me puede ayudar y últimamente se han estado utilizando R más y más, pero hay un análisis que hago algo de forma rutinaria que me ha dado problemas en R.

He de medidas repetidas de datos donde me gustaría ajustarse a la siguiente modelo: $$Delta = Day + Group + Day\times Group$$ where $Delta$ is the change from baseline, $Día$ is the number of days from the beginning of the study, and $$ es el grupo experimental. I ajuste de la varianza-covarianza de la matriz a cuenta para medidas repetidas (para este ejemplo estoy usando el compuesto de la simetría, pero la diferencia es el mismo utilizando otros varianza-covarianza de las matrices). Tengo los datos al final del post.

Si yo no incluye la interacción, puedo conseguir los análisis para que se ejecute como yo quiero que en el SAS y 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)

Obviamente, los resultados difieren en términos del denominador el DF, pero no tengo la intención de iniciar la discusión (a menos que la diferencia está causando el problema). Cuando agrego la interacción en SAS, todo está 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 de 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 quejan de que el ajuste está en singular, pero SAS no?

Aquí están algunos de los datos falsos que son representativos de los datos con los que trabajo (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í tienes 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;

3voto

Sam Dickson Puntos 451

Esta es sólo una respuesta parcial hasta el momento. Parece que el principal culpable es el hecho de que no todos los grupos tienen los datos en todos los puntos en el tiempo, que es lo que hace que la matriz a ser singular cuando se incluye el término de interacción, por lo que puede obtener resultados similares mediante la ejecución de:

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 acerca de qué SAS para obtener alrededor de la falta de información en determinados puntos en el tiempo que aún queda, pero habiendo descubierto el principal desencadenante me ayudará a avanzar hacia una resolución feliz.

0voto

RandomWhiteTrash Puntos 155

Recibí un mensaje de error similar y finalmente descubrí que tenía demasiados datos faltantes en la variable y de mi conjunto de datos para ciertas categorías de factores. Cuando reemplacé muchas de mis NA con números (afortunadamente en mi caso se suponía que esos NA eran ceros), el modelo funcionó bien.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by: