nlmer longitudinal data

Question

"nlme" paketindeki R Orthodont veri kümesiyle çalışıyorum. Bir göz atmak için sadece install.packages("nlme");library(nlme);head(Orthodont) kullanın. Veri seti, zaman içinde 27 çocukta ölçülen hipofiz ve pterjiyomaxiller fissür arasındaki mesafeden oluşmaktadır. lme4 paketinin kullanılması Doğrusal olmayan karma efekt modelini, fonksiyonel formum olarak bir lojistik eğri kullanarak takabilirim. Ben Asimptota olmasını seçebilir ve orta nokta olarak rastgele etkilere

cinsiyet bu parametreleri değişirse ne Gerçekten bilmek istiyorum

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) 

girdi. Ne yazık ki çevrimiçi örnekler hem konu hem de grup örneklerini içermez. Bu lme4 paketi ile mümkün mü?

Answer 1

Özel model formülasyonu ve eğimi için bir işlev oluşturarak böyle bir şey yapmanın mümkün olduğuna inanıyorum. Standart SSlogis fonksiyonu aşağıdaki biçimde bir lojistik fonksiyonunu kullanır:

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

yerine SSlogis çağırmak, kendi ihtiyaçlarını karşılamak için yukarıdaki ifadeyi değiştirebilir. Cinsiyetin sabit bir etki yaratıp yaratmadığını görmek istediğinize inanıyorum. İşte Asym2bir cinsiyet spesifik Asım alt popülasyon etkisi değiştirmek için örnek kod geçerli:

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

cinsiyet etkisinin tanıtılması Oldukça tipik yolu ikili kodlama olduğunu.

# 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"]),] 

Sonra özelleştirilmiş modeli sığabilecek: Ben Cinsiyet kodlu bir ikili -Değişken Seks dönüştürecek

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

Ne olur ki zaman Cinsiyet == 0 (Erkek), model değerlere ulaşır: h aslında standart SSlogis-fonksiyonu formudur. Ancak, ikili anahtarı şimdi olduğu takdirde Cinsiyet == 1 (Kadın):

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

yüzden yaş arttıkça elde asimptotik seviyesi aslında Asım + Asym2 değil, sadece Asım olduğuna Kadın bireyler için. Ayrıca, Asym2 için yeni bir rasgele efekt belirtmediğime dikkat edin. Asym, cinsiyete özgü olmadığı için, kadın bireyler bireylerinde Asym -term'den dolayı asimptotik düzeylerde varyansa sahip olabilirler. Model uyum: - 3,493 = 29.882:

> 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 

kadın hastalar 'yaş' arttıkça 'mesafeye' biraz daha düşük değerlere ulaşır görünüyor ki, bir cinsiyete özgü bir etkiye (t -2,859) olabilir gibi görünüyor 26.389

Sadece bunun iyi bir/en iyi model olduğunu aklınızdan çıkarmayacağım, sadece lme4 numaralı lineer olmayan modelleri özelleştirerek nasıl devam edebileceğinizi gösteriyorum.

# 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 

Grafik çıkışı:

sen (How do I extract lmer fixed effects by observation? doğrusal modeller için görselleştirme benzer şekilde) doğrusal olmayan sabit etkilerini ayıklamak istiyorsanız model için Görselleştirmeler müdahalesi biraz gerektirir

Yukarıdaki resimde, Dişi (F ##) bireylerin, Erkek (M ##) muadillerinden (siyah çizgiler) biraz daha düşük olduğuna dikkat edin. Örneğin. M10 < -> Orta alan panellerinde F10 farkı.

artığı ve belirtilen model bazı özelliklerini incelemek için rastgele bir etki. Bireysel M13 biraz zor görünüyor.