1
$\begingroup$

This is a follow-up question to Which model for my data? (testing the differences in slope for three groups).

The solution from there works (big thanks to Heteroskedastic Jim!), but I have a problem with a specific data set. Maybe someone can enlighten me why I get stuck.

Here is an example that works:

library(nlme)
library(emmeans)

Input = ("
Group   Time    Size
         A  1   1.08152
         A  2   1.10589
         A  3   1.13292
         B  1   1.04597
         B  2   1.05763
         B  3   1.07023
         B  4   1.08612
         B  5   1.10059
         B  6   1.11589
         B  7   1.13143
         B  8   1.14741
         B  9   1.16721
         B  10  1.18288
         C  1   1.04777
         C  2   1.06145
         C  3   1.07484
         C  4   1.08908
         C  5   1.10346
         C  6   1.11866
         C  7   1.13375
         C  8   1.14931
         C  9   1.16563
         C  10  1.18294
         ")
dat = read.table(textConnection(Input),header=TRUE)

This constructs the model:

(m1 <- gls(Size ~ Time * Group, dat, correlation = corAR1(form = ~ Time | Group), weights = varIdent(form = ~ 1 | I(Group == "A"))))

And this provides me with the p-values for slope differences:

pairs(emtrends(m1, ~ Group, var = "Time", df = Inf, options = get_emm_option("emmeans")))

Now the data set where I get stuck:

Input = ("
Group   Time    Size
         A  1   1.6210
         A  2   2.1118
         A  3   2.6026
         A  4   3.0934
         B  1   0.9162
         B  2   1.2122
         B  3   1.5082
         B  4   1.8042
         B  5   2.1002
         B  6   2.3962
         B  7   2.6922
         B  8   2.9882
         B  9   3.2842
         B  10  3.5802
         C  1   0.82701
         C  2   1.13441
         C  3   1.44181
         C  4   1.74921
         C  5   2.05661
         C  6   2.36401
         C  7   2.67141
         C  8   2.97881
         C  9   3.28621
         C  10  3.59361
         ")
dat = read.table(textConnection(Input),header=TRUE)

When I construct the above model with this specific data

(m1 <- gls(Size ~ Time * Group, dat, correlation = corAR1(form = ~ Time | Group), weights = varIdent(form = ~ 1 | I(Group == "A"))))

I get this error message:

Error in glsEstimate(object, control = control) : computed "gls" fit is singular, rank 6

I have tried analyzing the data in SPSS, but I also got stuck there.

So my question is: where is the problem with my data and what can I do to solve it?

Improve this question
$\endgroup$
4
  • $\begingroup$ So you kept the same model but changed the data, and the model cannot be estimated anymore? $\endgroup$
    – user2974951
    Commented Sep 19, 2019 at 13:11
  • $\begingroup$ Yes, this is correct. $\endgroup$
    – Kardashev3
    Commented Sep 20, 2019 at 7:38
  • $\begingroup$ Is time a factor or a variable? $\endgroup$
    – Sextus Empiricus
    Commented Sep 22, 2019 at 11:47
  • $\begingroup$ Time is a variable. $\endgroup$
    – Kardashev3
    Commented Sep 23, 2019 at 9:21

1 Answer 1

Reset to default
3
+50
$\begingroup$

If you plot Size against Time for each Group you will find that the points all lie on a straight line. Since you are fitting a model with the interaction between Time and Group you get a perfect fit overall which is what the software is telling you. Without knowing more about the process which generates your data it is impossible to say what implications this has for your scientific question.

Improve this answer
$\endgroup$
5
  • $\begingroup$ I think this doesn't explain it completely. If it was just a perfect fit, the fit could be computed and given out with variance zero, no reason for the software to report an error. I'd have thought that the error most likely means that the model is overparametrized, but this is also weird because then why does it work with the other dataset? There is a tiny difference in covariates though, an A 4 case is added, but I don't quite understand how this causes singularity (which is usually a feature of the model and design, not the response). $\endgroup$
    – Christian Hennig
    Commented Sep 24, 2019 at 18:37
  • $\begingroup$ @Lewian I fitted without the autocorrelation structure to keep it simple and then you get fitted results but a warning about perfect fit which does what you suggest it ought to. I suspect the added complexity of the gls() fit makes it fail. $\endgroup$
    – mdewey
    Commented Sep 25, 2019 at 8:51
  • $\begingroup$ What do you think, is there a way to test the three groups for differences in slope (and maybe intercept) using this or another suitable model (considering autocorrelation)? What baffles me is that R simply throws this error and stops. $\endgroup$
    – Kardashev3
    Commented Sep 25, 2019 at 9:01
  • 1
    $\begingroup$ You have a perfect fit without autocorrelation so adding it is not going to work. The model is not appropriate for this particular data-set. $\endgroup$
    – mdewey
    Commented Sep 25, 2019 at 10:10
  • $\begingroup$ @mdewey, yours is the correct answer. I dug into the table with the provided values and found the real problem. I was under the impression that I was given raw data, but instead the data provided resulted from a linear equation. Now it makes perfect sense why all the points are on a line :/. After digging to the raw data itself the model works... Mea culpa! Just out of curiosity: I am testing for differences in slope, could you give me a hint how to test for differences in intercept in the given model? $\endgroup$
    – Kardashev3
    Commented Sep 25, 2019 at 14:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.