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 fe `'combined'` — y = f + (a +bf) 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`	Количество итераций. Это может быть скаляром или трехэлементным вектором. Средства управления, сколько итераций выполняется для каждой из трех фаз алгоритма: симулированный отжиг полный размер шага уменьшаемый размер шага Значением по умолчанию является `[150 150 100]`. Скаляр распределяется на этих трех фазах в тех же пропорциях как значение по умолчанию.
`NMCMCIterations`	Количество итераций Цепи Маркова Монте-Карло (MCMC). Это может быть скаляром или трехэлементным вектором. Средства управления, сколько из трех различных типов обновлений MCMC выполняется во время каждой фазы основной итерации: полное многомерное обновление одно координатное обновление несколько координируют обновление Значением по умолчанию является `[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 = ADESIGNBETA + BDESIGNB 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 — Индивидуум взвесил остаточные значения

Примеры

свернуть все

Нелинейная модель Смешанных Эффектов со стохастическим алгоритмом EM

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

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

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}, Σ) = \int p (y | β, b, σ^{2}) p (b | Σ) d b$

где

y является данными об ответе
β вектор из коэффициентов населения
σ² остаточное отклонение
∑ ковариационная матрица для случайных эффектов
b является набором ненаблюдаемых случайных эффектов

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

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

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

Ссылки

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

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

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

nlinfit | nlpredci | nlmefit

Темы

Введен в R2010a

Документация