Exponenta Banner
exponenta event banner

nlmefitsa

Подбор нелинейной модели смешанных эффектов с алгоритмом stochastic 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 соответствует модели, где каждый параметр модели является суммой соответствующего фиксированного и случайного эффекта, а ковариационная матрица случайных эффектов является диагональной, то есть некоррелированными случайными эффектами.

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-на-ч из n наблюдений по переменным предсказателя h.

Y

Вектор ответов n-на-1.

GROUP

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

V

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

MODELFUN

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

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

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

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

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

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

  • VFUN - Либо

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

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

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

'FEParamsSelect'

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

'FEConstDesign'

Матрица проектирования 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'

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

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

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

'CovPattern'

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

Либо укажите PAT как вектор 1 на 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), или эквивалентно log (y) = log (f) + a * e

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

  • для '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' Требуется Toolbox™ оптимизации.

'Options'

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

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

  • 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 - Критерий информации Акайке (пуст, если 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 имеет значение false)

  • 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 'β, start2, Λ) =∫p (y' β, b, start2) p (b 'Λ) db

где

  • y - данные ответа

  • β - вектор коэффициентов популяции;

  • start2 - остаточная дисперсия

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

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

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

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

  1. Моделирование: Генерировать смоделированные значения случайных эффектов b из задней плотности p (b 'Λ), учитывая текущие оценки параметров.

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

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

Ссылки

[1] Делен, Б., М. Лавилль и Э. Мулайнс, «Сходимость стохастической аппроксимационной версии алгоритма ЕМ». Анналы статистики, 27, 94-128, 1999.

[2] Mentré, F. и M. Lavielle, «Стохастические алгоритмы ЕМ в популяционном анализе PKPD». Американская конференция по фармакометрии, 2008 год.

См. также

| |

Темы

Представлен в R2010a

Документация по инструментам для статистического и машинного обучения

Поддержка