Exponenta Banner

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

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 — 1 p вектором из параметров модели.

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

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

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

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

  • VFUN — Также

    • 1 g вектором из специфичных для группы предикторов для одной группы, соответствуя одной строке V

    • n-by-g матрица, где k-th строка VFUN V(i, :) если k-th наблюдение находится в группе 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 в виде набора начальных значений.

Аргументы name-value

По умолчанию, 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, где BDESIGNB случайные компоненты элементов p PHI. Эта матрица должна состоять из 0s и 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'. Каждая модель задает ошибку стандартная нормальная переменная (Gaussian) 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'. Использование 'fminunc' требует Optimization Toolbox™. 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уникальный(group)) + numelY)

  • sebeta — Стандартные погрешности для BETA (пустой, если ComputeStdErrors параметр имеет свое значение по умолчанию лжи),

  • 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 () функция на правой стороне является нормальной (Гауссовой) функцией правдоподобия, которая может зависеть от ковариантов.

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

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

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

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

Ссылки

[1] Delyon, B., М. Лэвилл и Э. Мулайнс, "Сходимость стохастической версии приближения алгоритма EM". Летопись Статистики, 27, 94-128, 1999.

[2] Mentré, F. и М. Лэвилл, "Стохастические Алгоритмы EM в Населении исследования PKPD". Американская Конференция по Фармакометрикам, 2008.

Смотрите также

| |

Темы

Введен в R2010a

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

Поддержка