R: как извлечь информацию из сводной модели соответствия

library(nlme)
fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
            data = Loblolly,
            fixed = Asym + R0 + lrc ~ 1,
            random = Asym ~ 1,
            start = c(Asym = 103, R0 = -8.5, lrc = -3.3))
> summary(fm1)
Nonlinear mixed-effects model fit by maximum likelihood
  Model: height ~ SSasymp(age, Asym, R0, lrc) 
 Data: Loblolly 
       AIC      BIC    logLik
  239.4856 251.6397 -114.7428

Random effects:
 Formula: Asym ~ 1 | Seed
            Asym  Residual
StdDev: 3.650642 0.7188625

Fixed effects: Asym + R0 + lrc ~ 1 
         Value Std.Error DF   t-value p-value
Asym 101.44960 2.4616951 68  41.21128       0
R0    -8.62733 0.3179505 68 -27.13420       0
lrc   -3.23375 0.0342702 68 -94.36052       0
 Correlation: 
    Asym   R0    
R0   0.704       
lrc -0.908 -0.827

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-2.23601930 -0.62380854  0.05917466  0.65727206  1.95794425 

Number of Observations: 84
Number of Groups: 14 

Я заинтересован в извлечении информации из итоговых результатов подгонки NLME.

Я хотел бы извлечь

  1. StdDev для случайных эффектов (т.е. StdDev для Asym, который = 3,65). Для этого я пробовал fm1$apVar, но не повезло.
  2. Оценки параметров фиксированных эффектов (например, Asym = 101,44960, R0 = -8,62733 и т. Д.), Которые можно извлечь с помощью fixef(fm1)
  3. Std. Ошибка фиксированных эффектов (т.е. 2.46, 0.317, 0.034). Для этого я пробовал sqrt(diag(fm1$varFix)), но эти значения не совсем соответствуют значениям в столбце Std.Error с фиксированными эффектами?
  4. logLikelihood (например, -114,7428, которое можно извлечь с помощью fm1$logLik)
  5. Остатки (например, 0,7188625, которые можно извлечь с помощью fm1$Residuals)

Моя конечная цель - подогнать несколько моделей и сохранить их соответствующие итоговые оценки в виде организованной data.frame.

fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
            data = Loblolly,
            fixed = Asym + R0 + lrc ~ 1,
            random = Asym ~ 1,
            start = c(Asym = 103, R0 = -8.5, lrc = -3.3))

fm2 <- nlme(height ~ SSasymp(age, Asym, R0, lrc),
            data = Loblolly,
            fixed = Asym + R0 + lrc ~ 1,
            random = Asym ~ 1,
            start = c(Asym = 103, R0 = -5.4, lrc = -3.3))

summary(fm1)
summary(fm2)

mylist = list(NULL, summary(fm1), NULL, summary(fm2), NULL, NULL)

Предположим, мой объект списка выглядит как mylist. Теперь я хотел бы создать data.frame, который выглядит так:

model    FixedAsym    FixedAsymStdError   FixedR0      ...     Residual
 1       101.44960        2.4616951       -8.62733            0.7188625
 2       101.44934        2.4616788       -8.62736     ...    0.7188625

И для создания этого data.frame (количество строк соответствует количеству сводок моделей, которые у меня есть в mylist) мне необходимо систематически извлекать эти значения (пронумерованные 1-5) из сводных результатов модели.


person Adrian    schedule 21.06.2017    source источник
comment
Создайте функцию, которая возвращает именованный вектор всей статистики, которую вы хотите вернуть для одной модели, а затем используйте sapply, чтобы запустить ее по списку моделей.   -  person lmo    schedule 21.06.2017
comment
@lmo Звучит хорошо. У вас есть идеи, как я могу извлечь файл std. ошибки фиксированных эффектов?   -  person Adrian    schedule 21.06.2017
comment
Вы можете использовать tidy и glance из пакета broom, чтобы получить большинство оценок модели, которые вам нужны. Альтернативой tidy является просто вытащить tTable из итогового вывода и извлечь оттуда то, что вы хотите (например, SE); например, summary(fm1)$tTable.   -  person aosmith    schedule 21.06.2017
comment
И логарифм правдоподобия, и остатки имеют функции извлечения: logLik & residuals. Вы можете сопоставить фиксированные эффекты ул. ошибки, возвращаемые sqrt(diag(vcov(fm1))) (или sqrt(diag(fm1$varFix), но опять же, лучше использовать функции извлечения) в итоговый вывод, если вы не корректируете их summary(fm1, adjustSigma = FALSE) или не соответствуете своей модели с помощью method="REML" - см. ?nlme:::summary.lme   -  person user20650    schedule 21.06.2017


Ответы (1)


Вот еще несколько штук ...

as.numeric(VarCorr(fm1)[,2])
# [1] 3.6506418 0.7188625

summary(fm1)$tTable[,2]
#       Asym         R0        lrc 
# 2.46169512 0.31795045 0.03427017 

# looks like you don't need this one anymore, but here's a way of getting it
summary(fm1)$corFixed
#            Asym         R0        lrc
# Asym  1.0000000  0.7039498 -0.9077793
# R0    0.7039498  1.0000000 -0.8271022
# lrc  -0.9077793 -0.8271022  1.0000000

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

person Matt Tyers    schedule 21.06.2017