dati longitudinali nlmer

Question

Ho lavorato con il set di dati R Orthodont nel pacchetto "nlme". Basta usare install.packages("nlme");library(nlme);head(Orthodont) per dare un'occhiata. Il set di dati comprende la distanza tra l'ipofisi e la fessura pterygomaxillary misurata in 27 bambini nel tempo. Utilizzo del pacchetto lme4 Posso adattare un modello di effetti misti non lineari utilizzando una curva logistica come forma funzionale. Posso scegliere di avere l'asintoto e punto medio entrato come effetti casuali

nm1 <- nlmer(distance ~ SSlogis(age,Asym, xmid, scal) ~ (Asym | Subject) + (xmid | Subject), Orthodont, start = c(Asym =25,xmid = 11, scal = 3), corr = FALSE,verb=1) 

Quello che vorrei davvero sapere è se il sesso cambia questi parametri. Sfortunatamente gli esempi online non includono sia l'oggetto che gli esempi di gruppo. Questo è possibile anche con il pacchetto lme4?

Answer 1

Credo sia possibile fare una cosa simile creando una funzione per la formulazione del modello personalizzata e il suo gradiente. La funzione standard SSlogis utilizza una funzione logistica della forma seguente:

f(input) = Asym/(1+exp((xmid-input)/scal)) # as in ?SSlogis 

Invece di chiamare i SSlogis, è possibile modificare la dichiarazione di cui sopra in base alle proprie esigenze. Credo che vorresti vedere se il genere ha effetto su un effetto fisso. Ecco codice di esempio per la modifica di un genere specifico Asym effetto sottopopolazione in Asym2:

# Just for loading the data, we will use lme4 for model fitting, not nlme 
library(nlme) 
library(lme4) 
# Careful when loading both nlme and lme4 as they have overlap, strange behaviour may occur 

# A more generalized form could be taken e.g. from http://en.wikipedia.org/wiki/Generalised_logistic_curve 
# A custom model structure: 
Model <- function(age, Asym, Asym2, xmid, scal, Gender) 
{ 
    # Taken from ?SSlogis, standard form: 
    #Asym/(1+exp((xmid-input)/scal)) 
    # Add gender-specific term to Asym2 
    (Asym+Asym2*Gender)/(1+exp((xmid-age)/scal)) 
    # Evaluation of above form is returned by this function 
} 

# Model gradient, notice that we include all 
# estimated fixed effects like 'Asym', 'Asym2', 'xmid' and 'scal' here, 
# but not covariates from the data: 'age' and 'Gender' 
ModelGradient <- deriv(
    body(Model)[[2]], 
    namevec = c("Asym", "Asym2", "xmid", "scal"), 
    function.arg=Model 
) 

Un modo piuttosto tipica di introdurre un effetto genere è con codifica binaria. Io trasformare la Sex -variable ad un codice binario Sesso:

# Binary coding for the gender 
Orthodont2 <- data.frame(Orthodont, Gender = as.numeric(Orthodont[,"Sex"])-1) 
#> table(Orthodont2[,"Gender"]) 
# 0 1 
#64 44 
# Ordering data based on factor levels so they don't mix up paneling in lattice later on 
Orthodont2 <- Orthodont2[order(Orthodont2[,"Subject"]),] 

posso quindi montare il modello personalizzato:

# Fit the non-linear mixed effects model 
fit <- nlmer(
    # Response 
    distance ~ 
    # Fixed effects 
    ModelGradient(age = age, Asym, Asym2, xmid, scal, Gender = Gender) ~ 
    # replaces: SSlogis(age,Asym, xmid, scal) ~ 
    # Random effects 
    (Asym | Subject) + (xmid | Subject), 
    # Data 
    data = Orthodont2, 
    start = c(Asym = 25, Asym2 = 15, xmid = 11, scal = 3)) 

Quello che succede è che quando Sesso == 0 (Maschio), il modello raggiunge i valori:

(Asym+Asym2*0)/(1+exp((xmid-age)/scal)) = (Asym)/(1+exp((xmid-age)/scal)) 

whic h è in realtà il modulo standard SSlogis. Tuttavia, v'è ora l'interruttore binario che se Sesso == 1 (femmina):

(Asym+Asym2)/(1+exp((xmid-age)/scal)) 

modo che il livello asintotico che otteniamo con l'aumentare dell'età è in realtà Asym + Asym2, non solo Asym, per individui di sesso femminile.

Si noti inoltre che non si specifica un nuovo effetto casuale per Asym2. Perché Asym non è specifico per il genere, le persone di sesso femminile possono anche avere varianza nei loro livelli asintotici di individui a causa di Asim -termine. Modello fit:

> summary(fit) 
Nonlinear mixed model fit by the Laplace approximation 
Formula: distance ~ ModelGradient(age = age, Asym, Asym2, xmid, scal,  Gender = Gender) ~ (Asym | Subject) + (xmid | Subject) 
    Data: Orthodont2 
    AIC BIC logLik deviance 
268.7 287.5 -127.4 254.7 
Random effects: 
Groups Name Variance Std.Dev. 
Subject Asym 7.0499 2.6552 
Subject xmid 4.4285 2.1044 
Residual  1.5354 1.2391 
Number of obs: 108, groups: Subject, 27 

Fixed effects: 
     Estimate Std. Error t value 
Asym 29.882  1.947 15.350 
Asym2 -3.493  1.222 -2.859 
xmid  1.240  1.068 1.161 
scal  5.532  1.782 3.104 

Correlation of Fixed Effects: 
     Asym Asym2 xmid 
Asym2 -0.471    
xmid -0.584 0.167  
scal 0.901 -0.239 -0.773 

Sembra che ci potrebbe essere un effetto specifico genere (t -2,859), in modo che i pazienti di sesso femminile sembrano raggiungere valori un po 'più bassi di 'distanza' con l'aumentare della 'età': 29,882-3,493 = 26.389

Non sto suggerendo necessariamente che questo è un modello buono/migliore, solo mostrando come si può andare avanti con la personalizzazione dei modelli non lineari in lme4. Effetti grafici per il modello richiedono un po 'di ritocchi se si desidera estrarre gli effetti non lineari fisse (in modo simile alla visualizzazione per i modelli lineari in How do I extract lmer fixed effects by observation?): uscita

# Extracting fixed effects components by calling the model function, a bit messy but it works 
# I like to do this for visualizing the model fit 
fixefmat <- matrix(rep(fixef(fit), times=dim(Orthodont2)[1]), ncol=length(fixef(fit)), byrow=TRUE) 
colnames(fixefmat) <- names(fixef(fit)) 
Orthtemp <- data.frame(fixefmat, Orthodont2) 
attach(Orthtemp) 
# see str(Orthtemp) 
# Evaluate the function for rows of the attached data.frame to extract fixed effects corresponding to observations 
fix = as.vector(as.formula(body(Model)[[2]])) 
detach(Orthtemp) 

nobs <- 4 # 4 observations per subject 
legend = list(text=list(c("y", "Xb + Zu", "Xb")), lines = list(col=c("blue", "red", "black"), pch=c(1,1,1), lwd=c(1,1,1), type=c("b","b","b"))) 
require(lattice) 
xyplot(
    distance ~ age | Subject, 
    data = Orthodont2, 
    panel = function(x, y, ...){ 
     panel.points(x, y, type='b', col='blue') 
     panel.points(x, fix[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='black') 
     panel.points(x, fitted(fit)[(1+nobs*(panel.number()-1)):(nobs*(panel.number()))], type='b', col='red') 
    }, 
    key = legend 
) 

# Residuals 
plot(Orthodont2[,"distance"], resid(fit), xlab="y", ylab="e") 

# Distribution of random effects 
par(mfrow=c(1,2)) 
hist(ranef(fit)[[1]][,1], xlab="Random 'Asym'", main="") 
hist(ranef(fit)[[1]][,2], xlab="Random 'xmid'", main="") 
# Random 'xmid' seems a bit skewed to the right and may violate normal distribution assumption 
# This is due to M13 having a bit abnormal growth curve (random effects): 
#   Asym  xmid 
#M13 3.07301310 3.9077583 

Grafica:

Nota come in figura sopra gli individui di sesso femminile (F ##) sono leggermente inferiori rispetto ai loro omologhi maschi (M ##) (linee nere). Per esempio. M10 < -> Differenza F10 nei pannelli dell'area centrale.

Residui ed effetti casuali per osservare alcune caratteristiche del modello specificato. L'individuo M13 sembra un po 'complicato.