0

Without providing a reproducible example (because that converged, and my data did not), my data is in a long format and includes 104 patients of which 304 measurements were taken. Most only have 2 measurements (n=68), followed by 5 (n=32), 4 (n=3), 3 (n=1), and 6 measurements (n=1). Some repeated measurements are taken fairly short after each other (for example for 74 of the 304 measurements the time in between was ~0.5 years), where for some patients with only two measurements the time in between is 5 years. The average follow-up time is about 3.5 years.

So far I've modeled a random intercept and slope model:

library(nlme)

ris <- 
  lme(fixed=outcome~ 1 + fu_time + age*fu_time + sex*fu_time +  smoking*fu_time + 
            obesity*fu_time + diab*fu_time + hypt*fu_time + hyperchol*fu_time + 
            ckd*fu_time, 
      random=~1 + fu_time|patid, 
      data=data, 
      na.action="na.omit",
      method="ML")

To further take into account the differing time in between subsequent measurements but keeping the increasing variances over time, I've tried to specify a model with (1) a random intercept, (2) a continuous AR1 correlation structure, and (3) heterogeneous variances:

ri_cAR1_hetvar <- 
  update(ris, 
         random=~1|patid, 
         correlation=corCAR1(form=~1|patid),
         weight=varIdent(form=~1|fu_time), 
         control=lmeControl(opt="optim")) # does not converge.

This model does not converge. As per advice of Ben Bolker previously here is some output from running debug(nlme:::logLik.reStruct) followed by the following commands:

### Running str(object)
List of 1
 $ patid: 'pdLogChol' num 0.0752
  ..- attr(*, "formula")=Class 'formula'  language ~1
  .. .. ..- attr(*, ".Environment")=<environment: 0x000000001ec31628> 
  ..- attr(*, "Dimnames")=List of 2
  .. ..$ : chr "(Intercept)"
  .. ..$ : chr "(Intercept)"
 - attr(*, "settings")= int [1:4] 0 1 0 4
 - attr(*, "class")= chr "reStruct"
 - attr(*, "plen")= Named int 1
  ..- attr(*, "names")= chr "patid"

### Running str(conLin)
List of 5
 $ Xy      : num [1:314, 1:22] 1 0.816 0.816 0.816 0.816 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:314] "1" "2" "3" "4" ...
  .. ..$ : chr [1:22] "(Intercept)" "(Intercept)" "fu_time" "age" ...
 $ dims    :List of 15
  ..$ N      : int 314
  ..$ ZXrows : int 314
  ..$ ZXcols : num 22
  ..$ Q      : int 1
  ..$ StrRows: num 125
  ..$ qvec   : Named num [1:3] 1 0 0
  .. ..- attr(*, "names")= chr [1:3] "patid" "" ""
  ..$ ngrps  : Named int [1:3] 104 1 1
  .. ..- attr(*, "names")= chr [1:3] "patid" "X" "y"
  ..$ DmOff  : Named num [1:3] 0 1 401
  .. ..- attr(*, "names")= chr [1:3] "" "patid" ""
  ..$ ncol   : Named num [1:3] 1 20 1
  .. ..- attr(*, "names")= chr [1:3] "patid" "" ""
  ..$ nrot   : Named num [1:3] 21 1 0
  .. ..- attr(*, "names")= chr [1:3] "" "" ""
  ..$ ZXoff  :List of 3
  .. ..$ patid: num [1:104] 0 6 8 13 15 18 23 28 30 32 ...
  .. ..$ X     : Named num 314
  .. .. ..- attr(*, "names")= chr "patid"
  .. ..$ y     : Named num 6594
  .. .. ..- attr(*, "names")= chr ""
  ..$ ZXlen  :List of 3
  .. ..$ patid: num [1:104] 6 2 5 2 3 5 5 2 2 2 ...
  .. ..$ X     : num 314
  .. ..$ y     : num 314
  ..$ SToff  :List of 3
  .. ..$ patid: num [1:104] 0 1 2 3 4 5 6 7 8 9 ...
  .. ..$ X     : Named num 229
  .. .. ..- attr(*, "names")= chr "patid"
  .. ..$ y     : Named num 2749
  .. .. ..- attr(*, "names")= chr ""
  ..$ DecOff :List of 3
  .. ..$ PXE_nr: num [1:104] 0 1 2 3 4 5 6 7 8 9 ...
  .. ..$ X     : Named num 125
  .. .. ..- attr(*, "names")= chr "patid"
  .. ..$ y     : Named num 2625
  .. .. ..- attr(*, "names")= chr ""
  ..$ DecLen :List of 3
  .. ..$ PXE_nr: num [1:104] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ X     : num 125
  .. ..$ y     : num 125
 $ logLik  :Class 'logLik' : 4.3 (df=178)
 $ sigma   : num 0
 $ auxSigma: num 0

### Running dput(object)
structure(list(patid= structure(-1.29387589179439, formula = ~1, Dimnames = list(
    "(Intercept)", "(Intercept)"), class = c("pdLogChol", "pdSymm", 
"pdMat"))), settings = c(0L, 1L, 0L, 4L), class = "reStruct", plen = c(patid= 1L))

I understand this question might be difficult to answer without a reproducible example, but posting anymore of my data is not really an option, and model works on the subsets of my data that I had prepared for this question. If anyone has any ideas how I could further diagnose this problem, the help would be very much appreciated!

Update I've noticed that the problem arises when I specify different variances per week by adding weight=varIdent(form=~1|fu_time), and nót per se by adding a continuous AR1 (that works fine). Might the problem be in how many parameters are estimated in the model with heterogenous variances? As in, I thought that that model still only just estimates phi, but now that I think of it it might estimate a parameter for every timepoint, which in my case is a continuous variable with a lot of levels...

4

0 回答 0