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

У меня есть многоуровневый набор данных повторных измерений около 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

r multi-level covariance nlme
person Annemarie    schedule 02.09.2016    source источник

Ответы (1)


arrow_upward
5
arrow_downward

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

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

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

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

подгонка против остатка

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

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

введите описание изображения здесь

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

Модель AR+random-slopes подходит нормально:

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"))

введите описание изображения здесь

person Ben Bolker    schedule 02.09.2016
comment
Большое спасибо за Вашу помощь. Да, данные такие колючие, и вы совершенно правы, логарифмическое преобразование делает их лучше. Я не могу понять, как добавить изображения в эти комментарии, но остатки всего набора данных в целом выглядят как ваши. Я изменил оптимизатор на nlminb, но теперь не могу добиться сходимости модели. Не могли бы вы дать еще какой-нибудь совет? Большое спасибо, Аннемари nlminb проблема, код ошибки сходимости = 1 сообщение = достигнут предел итерации без сходимости (10) - person Annemarie; 02.09.2016
comment
см. ?lmeControl (особенно maxIter, msMaxIter), хотя это может и не решить проблему. - person Ben Bolker; 02.09.2016
comment
Большое спасибо @Бен Болкер - person Annemarie; 02.09.2016
comment
Хотя это приветствуется, StackOverflow не рекомендует использовать комментарии, чтобы сказать спасибо ; если этот ответ был полезен, вы можете проголосовать за него (если у вас достаточно репутации), и в любом случае, если он удовлетворительно отвечает на ваш вопрос, вам предлагается нажать на галочку, чтобы принять его. - person Ben Bolker; 02.09.2016