nlmefitsa

Подгонка нелинейной модели смешанных эффектов с стохастическим алгоритмом EM

Свернуть все на странице

Синтаксис

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0)
[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name',Value)

Описание

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0) подходит для нелинейной регрессионной модели смешанных эффектов и возвращает оценки фиксированных эффектов в BETA. По умолчанию nlmefitsa подходит для модели, где каждый параметр модели является суммой соответствующего фиксированного и случайного эффекта, и ковариационная матрица случайных эффектов является диагональной, то есть некоррелированными случайными эффектами.

The BETA, PSIи другие значения, которые эти возвраты функции являются результатом случайного (Монте-Карло) симуляция, предназначенного для сходимости к максимальным оценкам вероятности параметров. Поскольку результаты случайны, желательно изучить график симуляции к результатам, чтобы убедиться, что симуляция сходилась. Также может быть полезно запустить функцию несколько раз, используя несколько начальных значений или использовать 'Replicates' параметр для выполнения нескольких симуляций.

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name',Value) принимает одну или несколько пар имя/значение параметра, разделенных запятыми. Задайте Name внутри одинарные кавычки.

Входные параметры

Определения:

В следующем списке аргументов применяются следующие определения переменных:

  • n - количество наблюдений

  • h - количество переменных предиктора

  • m - количество групп

  • g - количество групповых переменных предиктора

  • p - количество параметров

  • f - количество фиксированных эффектов

X

Матрица n -by h n наблюдений за переменными h.

Y

Вектор n -by-1 ответов.

GROUP

Сгруппированная переменная, указывающая, к какой из m групп принадлежит каждое наблюдение. GROUP может быть категориальной переменной, числовым вектором, символьной матрицей с строками для имен групп, строковыми массивами или массивом ячеек из векторов символов.

V

Матрица m -by g g предикторных переменных для каждой из m групп в данных. Это значения предиктора, которые берут одно и то же значение для всех наблюдений в группе. Строки V упорядочены согласно GRP2IDX(GROUP). Используйте m -by g массив ячеек для V если любое из значений предиктора для группы варьируется в размере по группам. Задайте [] для V если предикторов групп нет.

MODELFUN

Указатель на функцию, которая принимает значения предиктора и параметры модели и возвращает подобранные значения. MODELFUN имеет форму YFIT = MODELFUN(PHI,XFUN,VFUN) с входными параметрами

  • PHI - вектор параметров модели p 1 байт .

  • XFUN - l -by - h массив переменных-предикторов, где

    • l равно 1, если XFUN является одной строкой X

    • l ni, если XFUN содержит строки X для одной группы размеров ni

    • l n, если XFUN содержит все строки X.

  • VFUN - Либо

    • Вектор g 1 байт предикторов для одной группы, соответствующий одной строке V

    • Матрица n -by g, где k -я строка VFUN является V(i,:) если k -е наблюдение находится в групповом i.

    Если V пуст, nlmefitsa вызывает MODELFUN с двумя входами.

MODELFUN возвращает вектор l -by-1 подобранных значений YFIT. Когда либо PHI или VFUN содержит одну строку, которая соответствует всем строкам в двух других входных параметрах. Для улучшенной производительности используйте 'Vectorization' пара имя/значение параметра (описано ниже), если MODELFUN может вычислять YFIT для более чем одного вектора параметров модели в одном вызове.

BETA0

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

Аргументы в виде пар имя-значение

По умолчанию nlmefitsa подходит для модели, где каждый параметр модели является суммой соответствующего фиксированного и случайного эффекта. Используйте следующие пары имя/значение параметра, чтобы подогнать модель с другим количеством или зависимостью от фиксированных или случайных эффектов. Используйте самое большее одно имя параметра с 'FE' префикс и одно имя параметра со 'RE' префикс. Обратите внимание, что некоторые варианты меняют способ nlmefitsa вызывает MODELFUN, как описано далее ниже.

'FEParamsSelect'

Вектор, определяющий, какие элементы вектора параметра модели PHI включите фиксированный эффект, как числовой вектор с элементами в 1: p или как 1-байтовый p логический вектор. Модель будет включать f фиксированных эффектов, где f является заданным количеством элементов.

'FEConstDesign'

A p -by f матрица проекта ADESIGN, где ADESIGN* BETA являются фиксированными компонентами p элементов PHI.

'FEGroupDesign'

Массив p -by f -by m, задающий другую матрицу проекта p -by f фиксированных эффектов для каждой из групп m.

'REParamsSelect'

Вектор, определяющий, какие элементы вектора параметра модели PHI включить случайный эффект, как числовой вектор с элементами в 1: p, или как 1-байтовый p логический вектор. Модель будет включать r случайных эффектов, где r - заданное количество элементов.

'REConstDesign'

A p -by r матрица проекта BDESIGN, где BDESIGN* B являются случайными компонентами p элементов PHI. Эта матрица должна состоять из 0s и 1s, с самое большее одним 1 в строку.

Модель по умолчанию эквивалентна установке обоих FEConstDesign и REConstDesign на eye(p), или установить оба FEParamsSelect и REParamsSelect в 1: p.

Дополнительные необязательные пары имя/значение параметров управляют итеративным алгоритмом, используемым для максимизации вероятности:

'CovPattern'

Задает r -by r логическую или числовую матрицу PAT который задает шаблон ковариационной матрицы случайных эффектов PSI. nlmefitsa вычисляет оценки для отклонений по диагонали PSI а также ковариации, которые соответствуют ненулевым в off-диагонали PAT. nlmefitsa ограничивает оставшиеся ковариации, т.е. те, которые соответствуют недиагональным нулям в PAT, чтобы быть нулем. PAT должна быть сочетанием в строке-столбце блочной диагональной матрицы и nlmefitsa добавляет ненулевые элементы в PAT по мере необходимости для создания такого шаблона. Значение по умолчанию PAT является eye(r), соответствующий некоррелированным случайным эффектам.

Кроме того, задайте PAT как вектор 1-by- r, содержащий значения в 1: r. В этом случае элементы PAT с равными значениями задают группы случайных эффектов, nlmefitsa оценивает ковариации только в группах и ограничивает ковариации между группами равными нулю.

'Cov0'

Начальное значение для ковариационной матрицы PSI. Должна быть r -by r положительно определенной матрицей. Если он пуст, значение по умолчанию зависит от значений BETA0.

'ComputeStdErrors'

true вычислить стандартные ошибки для оценок коэффициентов и сохранить их в выход STATS структура, или false (по умолчанию), чтобы опустить эти расчеты.

'ErrorModel'

Вектор символов или строковый скаляр, задающий форму термина ошибки. По умолчанию это 'constant'. Каждая модель определяет ошибку, используя стандартную нормальную (Гауссову) переменную e, значение функции f и один или два параметра a и b. Варианты:

  • 'constant' — y = f + a * e

  • 'proportional' — y = f + b * f * e

  • 'combined' — y = f + (a + b * f) * e

  • 'exponential' - y = f * exp (a * e), или эквивалентно логгирование (y) = журналу (f) + a * e

Если задан этот параметр, выход STATS.errorparam поле имеет значение

  • a для 'constant' и 'exponential'

  • b для 'proportional'

  • [a b] для 'combined'

'ErrorParameters'

Скалярный или двухэлементный вектор, задающий начальные значения для параметров модели ошибки. Это задает значения a, b или [a b] в зависимости от ErrorModel параметр.

'LogLikMethod'

Задает метод аппроксимации логарифмической правдоподобности. Варианты:

  • 'is' - Отбор проб важности

  • 'gq' - Гауссова квадратура

  • 'lin' - Линеаризация

  • 'none' - Опустить логарифмическое приближение (по умолчанию)

'NBurnIn'

Количество начальных итераций горения, во время которых оценки параметра не пересчитываются. Значение по умолчанию является 5.

'NChains'

Количество c цепей». Значение по умолчанию является 1. Установка c > 1 приводит c вычислению векторов моделируемых коэффициентов для каждой группы во время каждой итерации. Значение по умолчанию зависит от данных и выбирается таким образом, чтобы обеспечить около 100 групп по всем цепям.

'NIterations'

Количество итераций. Это может быть скаляр или трехэлементный вектор. Управляет, сколько итераций выполняется для каждой из трех фаз алгоритма:

  1. моделируемый отжиг

  2. полный размер шага

  3. уменьшенный размер шага

По умолчанию это [150 150 100]. Скаляр распределен между тремя фазами в тех же пропорциях, что и по умолчанию.

'NMCMCIterations'

Количество итераций марковской цепи Монте-Карло (MCMC). Это может быть скаляр или трехэлементный вектор. Управляет тем, сколько из трех различных типов обновлений MCMC выполняется в течение каждой фазы основной итерации:

  1. полное многомерное обновление

  2. обновление одной координаты

  3. обновление нескольких координат

По умолчанию это [2 2 2]. Скалярное значение обрабатывается как трехэлементный вектор со всеми элементами, равными скаляру.

'OptimFun'

Либо 'fminsearch' или 'fminunc', определение оптимизационной функции, которая будет использоваться в процессе оценки. По умолчанию это 'fminsearch'. Использование 'fminunc' требует Optimization Toolbox™.

'Options'

Структура, созданная вызовом statset. nlmefitsa использует следующие statset параметры:

  • 'DerivStep' - Относительное различие, используемая в вычислении градиента конечного различия. Может быть скаляром или вектором, чья длина является количеством параметров модели p. Значение по умолчанию является eps^(1/3).

  • Display - Level of display во время оценки.

    • 'off' (по умолчанию) - Не отображает никакой информации

    • 'final' - Отображает информацию после окончательной итерации алгоритма оценки

    • 'iter' - Отображает информацию при каждой итерации

  • FunValCheck

    • 'on' (по умолчанию) - Проверяйте на недопустимые значения (такие как NaN или Inf) из MODELFUN

    • 'off' - Пропустить эту проверку

  • OutputFcn - Указатель на функцию, заданный с помощью @, массив ячеек с указателями на функцию или пустой массив. nlmefitsa вызывает все выходные функции после каждой итерации. См. nlmefitoutputfcn.m (выходная функция по умолчанию для nlmefitsa) для примера выходной функции.

'ParamTransform'

Вектор p -значений, задающих функцию преобразования f() для каждого из p параметров: 

XB = ADESIGN*BETA + BDESIGN*B 
PHI = f(XB)
Каждый элемент вектора должен быть одним из следующих целочисленных кодов, определяющих преобразование для соответствующего значения PHI:

  • 0: PHI = XB (по умолчанию для всех параметров)

  • 1: log(PHI) = XB

  • 2: probit(PHI) = XB

  • 3: logit(PHI) = XB

'Replicates'

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

'Vectorization'

Определяет возможные размеры PHI, XFUN, и VFUN входные параметры в MODELFUN. Возможные значения:

  • 'SinglePhi' — MODELFUN является функцией (такой как решатель ОДУ), которая может только вычислять YFIT для одного набора параметров модели за раз, т.е. PHI должен быть одна строка вектора при каждом вызове. nlmefitsa вызывает MODELFUN в цикле, если необходимо использование одной PHI вектор и с XFUN содержащие строки для одного наблюдения или группы за раз. VFUN может быть одна строка, которая применяется ко всем строкам XFUN, или матрица со строками, соответствующими строкам в XFUN.

  • 'SingleGroup' — MODELFUN может принимать только входы, соответствующие одной группе в данных, т.е. XFUN должны содержать строки X из одной группы в каждом вызове. В зависимости от модели, PHI - одна строка, которая применяется ко всей группе или матрица с одной строкой для каждого наблюдения. VFUN - одна строка.

  • 'Full' — MODELFUN может принимать входы для нескольких векторов параметров и нескольких групп в данных. Либо PHI или VFUN может быть одна строка, которая применяется ко всем строкам XFUN, или матрица со строками, соответствующими строкам в XFUN. Использование этой опции может улучшить эффективность путем уменьшения количества вызовов до MODELFUN, но может потребовать MODELFUN для выполнения расширения синглтона на PHI или V.

Значение по умолчанию для 'Vectorization' является 'SinglePhi'. Во всех случаях, если V пуст, nlmefitsa вызывает MODELFUN с двумя входами.

Выходные аргументы

BETA

Оценки фиксированных эффектов

PSI

r -by r оцененная ковариационная матрица для случайных эффектов. По умолчанию r равно количеству p параметров модели.

STATS

Структура со следующими полями:

  • logl - максимальная логарифмическая правдоподобность для установленной модели; пустой, если LogLikMethod параметр имеет значение по умолчанию 'none'

  • rmse - Квадратный корень предполагаемого отклонения ошибок (вычисленный по шкале журнала для exponential модель ошибки)

  • errorparam - Предполагаемые параметры модели отклонения ошибок

  • aic - Информационный критерий Akaike (пуст, если logl пуст), рассчитывается как aic = –2 * logl + 2 * numParam, где

    • logl - максимальная логарифмическая правдоподобность.

    • numParam - количество подходящих параметров, включая степень свободы для ковариационной матрицы случайных эффектов, количество фиксированных эффектов и количество параметров модели ошибки.

  • bic - Байесовский информационный критерий (пуст, если logl пуст), рассчитывается как bic =-2* logl + журнал (M) * numParam

    • M количество групп.

    • logl и numParam определяются как в aic.

    Обратите внимание, что в некоторой литературе предполагается, что расчет bic должен быть, bic =-2* logl + журнал (N) * numParam, где N количество наблюдений. Чтобы настроить значение выхода, можно переопределить bic следующим образом: bic = bic - numel(unique(group)) + numel(Y)

  • sebeta - Стандартные ошибки для BETA (пустые, если ComputeStdErrors параметр имеет значение по умолчанию false)

  • covb - Оценочная ковариация оценок параметров (пустая, если ComputeStdErrors является ложным)

  • dfe - ошибка степеней свободы

  • pres - остатки населения (y-y_population), где y_population является население предсказанными значениями

  • ires - остатки населения (y-y_population), где y_population является индивидуальными предсказанными значениями

  • pwres - население взвешенных невязок

  • cwres - Условные взвешенные невязки

  • iwres - Отдельные взвешенные невязки

Примеры

свернуть все

Открыть скрипт

Загрузите выборочные данные.

load indomethacin

Подгонка модели для данных о концентрациях препарата индометацина в кровотоке шести субъектов в течение восьми часов.

model = @(phi,t)(phi(:,1).*exp(-phi(:,2).*t)+phi(:,3).*exp(-phi(:,4).*t));
phi0 = [1 1 1 1];
xform = [0 1 0 1]; % log transform for 2nd and 4th parameters
[beta,PSI,stats,br] = nlmefitsa(time,concentration,...
   subject,[],model,phi0,'ParamTransform',xform)

beta =

    0.8563
   -0.7950
    2.7744
    1.0772


PSI =

    0.0529         0         0         0
         0    0.0220         0         0
         0         0    0.4762         0
         0         0         0    0.0120


stats = 

  struct with fields:

          logl: []
           aic: []
           bic: []
        sebeta: []
           dfe: 57
          covb: []
    errorparam: 0.0809
          rmse: 0.0775
          ires: [66x1 double]
          pres: [66x1 double]
         iwres: [66x1 double]
         pwres: [66x1 double]
         cwres: [66x1 double]


br =

   -0.2255    0.0063    0.1600    0.1773   -0.3269    0.1157
    0.0350   -0.1384    0.0058    0.0431    0.0093   -0.0453
   -0.7557   -0.0550    0.8736   -0.7875    0.5304    0.1727
   -0.0010   -0.0198    0.0137   -0.0757    0.0478   -0.0076

Постройте график данных вместе с общим населением подгонки

clf
phi = [beta(1), exp(beta(2)), beta(3), exp(beta(4))];
h = gscatter(time,concentration,subject);
xlabel('Time (hours)')
ylabel('Concentration (mcg/ml)')
title('{\bf Indomethacin Elimination}')
xx = linspace(0,8);
line(xx,model(phi,xx),'linewidth',2,'color','k')

Постройте график отдельных кривых на основе оценок случайного эффекта.

for j=1:6
    phir = [beta(1)+br(1,j), exp(beta(2)+br(2,j)), ...
            beta(3)+br(3,j), exp(beta(4)+br(4,j))];
    line(xx,model(phir,xx),'color',get(h(j),'color'))
end

Алгоритмы

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

p(y|β,σ2,Σ)=∫p(y|β,b,σ2)p(b|Σ)db

где

  • y - данные отклика

  • β - вектор населений

  • σ2 - остаточное отклонение

  • ∑ является ковариационной матрицей для случайных эффектов

  • b - это набор ненаблюдаемых случайных эффектов

Каждая p () функция на правой стороне является нормальной (Гауссовой) функцией правдоподобия, которая может зависеть от ковариат .

Поскольку интеграл не имеет замкнутой формы, трудно найти параметры, которые максимизируют его. Делион, Лавиэль и Мулин [1] предложили найти максимальные оценки правдоподобия с помощью алгоритма максимизации ожиданий (EM), в котором шаг E заменяется стохастической процедурой. Они назвали свой алгоритм SAEM для Стохастического приближения EM. Они продемонстрировали, что этот алгоритм имеет желательные теоретические свойства, включая сходимость в практических условиях и сходимость к локальному максимуму функции правдоподобия. Их предложение предусматривает три этапа:

  1. Симуляция: Сгенерируйте моделируемые значения случайных эффектов, b из апостериорной p плотности (b |

  2. Стохастическое приближение: Обновите ожидаемое значение функции логарифмической правдоподобности, взяв ее значение от предыдущего шага и двигаясь частично к среднему значению логарифмической правдоподобности, рассчитанному из моделируемых случайных эффектов.

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

Ссылки

[1] Delyon, B., M. Lavielle, and E. Moulines, «Сходимость стохастической версии приближения алгоритма EM». Анналы статистики, 27, 94 - 128, 1999.

[2] Mentré, F., and M. Lavielle, «Stochastic EM Algorithms in Population PKPD analyses». Американская конференция по фармакометрии, 2008.

См. также

| |

Темы

Введенный в R2010a

Statistics and Machine Learning Toolbox документация

Поддержка

Для просмотра документации необходимо авторизоваться на сайте