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

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 как набор начальных значений.

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

По умолчанию 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. Эта матрица должна состоять из 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'

Задает метод для приближения loglikelihood. Выбор:

  • 'is' — Выборка важности

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

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

  • 'none' Не используйте loglikelihood приближение (значение по умолчанию)

'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).

  • Отображение Уровень отображения во время оценки.

    • 'off' (значение по умолчанию) — Отображения никакая информация

    • 'final' — Информация об отображениях после итоговой итерации алгоритма оценки

    • 'iter' — Информация об отображениях в каждой итерации

  • FunValCheck

    • 'on' (sdefault) — Проверяйте на недопустимые значения (такие как 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 является одной строкой.

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

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

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

BETA

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

PSI

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

STATS

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

  • logl — Максимизируемый loglikelihood для подобранной модели; пустой, если параметр LogLikMethod имеет свое значение по умолчанию 'none'

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

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

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

    • logl является максимизируемым loglikelihood.

    • 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 имеет свое значение по умолчанию лжи),

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

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

Ссылки

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

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

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

| |

Темы

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

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