структура ковариации для многоуровневого моделирования

Question

У меня есть многоуровневый набор данных повторных измерений около 300 пациентов, каждый из которых имеет до 10 повторных мер, предсказывающих рост тропонина. В наборе данных есть и другие переменные, но я не включил их здесь. Я пытаюсь использовать nlme для создания случайной кривой, случайной модели перехвата, где эффекты варьируются между пациентами, а эффект времени у разных пациентов различен. Когда я пытаюсь ввести структуру ковариации первого порядка, чтобы обеспечить корреляцию измерений из-за времени, я получаю следующее сообщение об ошибке.

Error in `coef<-.corARMA`(`*tmp*`, value = value[parMap[, i]]) : Coefficient matrix not invertible 

Я включил мой код и образец набора данных, и я был бы очень благодарен за любые слова мудрости.

#baseline model includes only the intercept. Random slopes - intercept varies across patients 
randomintercept <- lme(troponin ~ 1, 
         data = df, random = ~1|record_id, method = "ML", 
         na.action = na.exclude, 
         control = list(opt="optim")) 

#random intercept and time as fixed effect 
timeri <- update(randomintercept,.~. + day) 
#random slopes and intercept: effect of time is different in different people 
timers <- update(timeri, random = ~ day|record_id) 
#model covariance structure. corAR1() first order autoregressive covariance structure, timepoints equally spaced 
armodel <- update(timers, correlation = corAR1(0, form = ~day|record_id)) 
Error in `coef<-.corARMA`(`*tmp*`, value = value[parMap[, i]]) : Coefficient matrix not invertible 

данные:

record_id day troponin 
1 1 32 
2 0  NA 
2 1 NA 
2 2 NA 
2 3 8 
2 4 6 
2 5 7 
2 6 7 
2 7 7 
2 8 NA 
2 9 9 
3 0 14 
3 1 1167 
3 2 1935 
4 0 19 
4 1 16 
4 2 29 
5 0 NA 
5 1 17 
5 2 47 
5 3 684 
6 0 46 
6 1 45440 
6 2 47085 
7 0 48 
7 1 87 
7 2 44 
7 3 20 
7 4 15 
7 5 11 
7 6 10 
7 7 11 
7 8 197 
8 0 28 
8 1 31 
9 0 NA 
9 1 204 
10 0 NA 
10 1 19 

Answer 1

Вы можете подходите, если вы измените свой оптимизатор «nlminb» (или, по крайней мере, он работает с уменьшенным набором данных, вы писали).

armodel <- update(timers, 
       correlation = corAR1(0, form = ~day|record_id), 
       control=list(opt="nlminb")) 

Однако, если вы посмотрите на подобранную модели, вы увидите, у вас есть проблемы - оценочный параметр П1 является -1, а случайное перехватывать и наклон условие коррелирует с г = 0,998.

Я думаю, что проблема связана с природой данных. Большинство данных, по-видимому, находятся в диапазоне 10-50, но есть экскурсии на один или два порядка (например, 6 человек, до 45000). Может быть, трудно подгонять модель к данным. Я бы сильно предлагаю log-преобразование ваших данных; стандартный диагностический участок (plot(randomintercept)) выглядит следующим образом:

, тогда как фитинг на логарифмической шкале

rlog <- update(randomintercept,log10(troponin) ~ .) 
plot(rlog) 

несколько более разумным, хотя есть еще некоторые доказательства гетероскедастичности.

модели AR + случайные наклоны подходит OK:

ar.rlog <- update(rlog, 
       random = ~day|record_id, 
       correlation = corAR1(0, form = ~day|record_id)) 
## Linear mixed-effects model fit by maximum likelihood 
## ... 
## Random effects: 
## Formula: ~day | record_id 
## Structure: General positive-definite, Log-Cholesky parametrization 
##    StdDev Corr 
## (Intercept) 0.1772409 (Intr) 
## day   0.6045765 0.992 
## Residual 0.4771523  
## 
## Correlation Structure: ARMA(1,0) 
## Formula: ~day | record_id 
## Parameter estimate(s): 
##  Phi1 
## 0.09181557 
## ... 

Быстрый взгляд на intervals(ar.rlog) показывает, что доверительные интервалы по авторегрессии параметра (-0.52,0.65), поэтому он не может быть стоит оставить ...

с случайных склонов в модели гетероскедастичности уже не кажется проблематичным ...

plot(rlog,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth"))