structure de covariance pour la modélisation multiniveau

Question

Je dispose d'un ensemble de données multiniveau à mesures répétées d'environ 300 patients, chacune contenant jusqu'à 10 mesures répétées permettant de prédire l'élévation de la troponine. Il y a d'autres variables dans l'ensemble de données, mais je ne les ai pas incluses ici. J'essaie d'utiliser nlme pour créer une pente aléatoire, modèle d'interception aléatoire où les effets varient entre les patients, et l'effet du temps est différent chez les différents patients. Lorsque j'essaie d'introduire une structure de covariance de premier ordre pour permettre la corrélation des mesures dues au temps, j'obtiens le message d'erreur suivant. J'ai inclus mon code et un échantillon de l'ensemble de données, et je serais très reconnaissant pour tous les mots de sagesse.

#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 

données:

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

Vous pouvez adapter si vous changez optimiseur à « nlminb » (ou au moins cela fonctionne avec les données réduites vous définissez au courant).

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

Cependant, si vous regardez le modèle ajusté, vous verrez que vous avez des problèmes - le paramètre AR1 est estimé -1 et les termes d'interception au hasard et la pente sont en corrélation avec r = 0,998.

Je pense que le problème est lié à la nature des données. La plupart des données semblent se situer dans l'intervalle de 10 à 50, mais il y a des excursions d'un ou deux ordres de grandeur (par exemple, individu 6, jusqu'à environ 45 000). Il pourrait être difficile d'adapter un modèle à des données aussi épineuses. Je voudrais fortement suggérer log-transformer vos données; le tracé de diagnostic standard (plot(randomintercept)) ressemble à ceci:

alors que le montage sur l'échelle logarithmique

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

est un peu plus raisonnable, bien qu'il y ait encore des preuves de hétéroscédasticité.

AR + modèle aléatoire de ski OK convient:

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 
## ... 

Un rapide coup d'œil à intervals(ar.rlog) montre que les intervalles de confiance sur le paramètre autorégressif sont (-0.52,0.65), donc il ne peut pas être utile de garder ...

Avec les pentes aléatoires dans le modèle du hétéroscédasticité ne semble plus problématique ...

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