0

I have fitted an lme model using a factor variable studytreat. It appears in the model in two places, as a fixed effect, and also to model heterogeneity in the residual variance so that each level of studytreat has its own variance. This is the model:

modelLME <- lme(logkp ~ studytreat, 
                random = ~ -1 + STUDY | ptno, 
                weights = varIdent(form = ~1 | studytreat))

I have a problem which doesn't fill me with confidence in lme! If I label the levels in studytreat as 11, 12, 21, 24 (in that order) then the function works fine. However, if I switch the numbers around to get 11, 21, 12, 42 then when I call summary(modelLME)$apVar I get the error message:

[1] "Non-positive definite approximate variance-covariance"

Does anyone have any idea why this should make a difference?

Tom Morris
  • 1
  • 1
  • 2
    Please provide the original factor variable and show the code you used to change the levels. – Sven Hohenstein Dec 20 '13 at 17:23
  • This post should also be rephrased as a question, along the lines of asking why it makes a difference. – Ganesh Sittampalam Dec 20 '13 at 17:42
  • Does the model that "works" really work, i.e. does it provide a solution with reasonable confidence intervals on the variances (what is `summary(modelLME)$apVar`) ? I can imagine that a just-barely-stable model could get tipped over the line to a not-quite-stable model by switching the structure ... – Ben Bolker Dec 20 '13 at 17:54
  • Ganesh: thanks, I have done that. – Tom Morris Dec 23 '13 at 09:29
  • Sven: I create studytreat from two numerical variables 'STUDY' = 1 or 2 and 'pktrt' = 1,2 or 4, by studytreat = 10*STUDY + pktrt, then studytreat <- factor(studytreat). For the reversal the roles of 'STUDY' and 'pktrt' are switched. – Tom Morris Dec 23 '13 at 09:31
  • Ben: I'm not sure what I should be looking for in summary(modelLME)$apVar to see if the model is unstable. I think it is quite close to being overparameterised so you might be right. However, I wouldn't have thought a relabelling should have any effect. I have noted since my original post that the order of the levels changes in the relabelling, which I suppose might count as a change in "structure". Thanks for your help. – Tom Morris Dec 23 '13 at 09:37
  • @TomMorris Thanks for the information. Why did you change the factor? – Sven Hohenstein Dec 23 '13 at 12:05
  • Sven: When I wrote down the model I denoted the outcome y_ijk where j is treatment and k is study. I wanted the factor labels to be in the same order as the subscripts, i.e. treatment before study. It doesn't matter really but I thought the strange behaviour of nlme worth mentioning. – Tom Morris Dec 23 '13 at 12:15

0 Answers0