nlmer dữ liệu theo chiều dọc

Question

Tôi đã làm việc với bộ dữ liệu R Orthodont trong gói "nlme". Chỉ cần sử dụng install.packages("nlme");library(nlme);head(Orthodont) để xem. Tập dữ liệu bao gồm khoảng cách giữa tuyến yên và vết nứt pterygomaxillary đo được ở 27 trẻ em theo thời gian. Sử dụng gói lme4 tôi có thể vừa với mô hình hiệu ứng hỗn hợp phi tuyến bằng cách sử dụng đường cong hậu cần làm biểu mẫu chức năng của tôi. Tôi có thể chọn để có asymptote và trung điểm được nhập như là các hiệu ứng ngẫu nhiên

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) 

Điều tôi thực sự muốn biết là nếu giới tính thay đổi các tham số này. Thật không may, ví dụ trực tuyến không bao gồm cả ví dụ về chủ đề và nhóm. Điều này thậm chí có thể với gói lme4?

Answer 1

Tôi tin rằng có thể thực hiện một việc như vậy bằng cách tạo một hàm cho công thức xây dựng mô hình tùy chỉnh và độ dốc của nó. Các tiêu chuẩn SSlogis chức năng sử dụng một hàm logistic có dạng sau:

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

Thay vì gọi SSlogis, bạn có thể sửa đổi tuyên bố trên cho phù hợp với nhu cầu của bạn. Tôi tin rằng bạn sẽ muốn xem liệu giới tính có ảnh hưởng đến hiệu ứng cố định hay không. Đây là mã ví dụ cho sửa đổi một giới tính cụ thể Asym quần thể có hiệu lực vào 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 
) 

Một cách khá điển hình của việc giới thiệu một tác dụng giới là với mã nhị phân. Tôi sẽ làm thay đổi Sex -variable đến một nhị phân mã Giới tính:

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

tôi sau đó có thể phù hợp với những mô hình tùy chỉnh:

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

Điều gì xảy ra khi == Giới tính 0 (Nam), mô hình đạt được các giá trị:

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

whic h thực sự là dạng hàm SSlogis tiêu chuẩn. Tuy nhiên, hiện nay là việc chuyển đổi nhị phân rằng nếu == Giới tính 1 (Nữ):

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

để mức tiệm cận mà chúng ta đạt được khi độ tuổi tăng thực sự là Asym + Asym2, không chỉ Asym, dành cho cá nhân nữ.

Cũng lưu ý rằng tôi không chỉ định hiệu ứng ngẫu nhiên mới cho Asym2. Bởi vì Asym không dành riêng cho giới tính, cá nhân nữ cũng có thể có sự khác biệt ở mức độ tiệm cận cá nhân do Asym -term. Mô hình phù hợp:

> 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 

Hình như có thể có một tác động cụ thể giới tính (t -2,859), do đó bệnh nhân nữ dường như đạt giá trị thấp hơn một chút 'khoảng cách' là 'tuổi' tăng: 29,882-3,493 = 26.389

Tôi không nhất thiết cho rằng đây là mô hình tốt/tốt nhất, chỉ hiển thị cách bạn có thể tiếp tục với tùy chỉnh mô hình phi tuyến tính trong lme4. Visualizations cho mô hình đòi hỏi một chút mày mò nếu bạn muốn trích xuất các hiệu ứng phi tuyến tính cố định (trong thời trang tương tự như trực quan cho các mô hình tuyến tính trong How do I extract lmer fixed effects by observation?):

# 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 

Graphics đầu ra:

Lưu ý cách thức ở trên hình các cá thể Nữ (F ##) hơi thấp hơn so với các đối tác Nam (M ##) của chúng (đường màu đen). Ví dụ. M10 < -> Sự khác biệt F10 trong bảng điều khiển ở giữa.

Dư và các hiệu ứng ngẫu nhiên để quan sát một số đặc điểm của mô hình cụ thể. M13 cá nhân có vẻ hơi phức tạp.