Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

nlmefitsa

Ajuste el modelo de efectos mixtos no lineales con el algoritmo estocástico EM

Sintaxis

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

Descripción

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0) se ajusta a un modelo de regresión de efectos mixtos no lineales y devuelve estimaciones de los efectos fijos en.BETA De forma predeterminada, se ajusta a un modelo donde cada parámetro de modelo es la suma de un efecto fijo y aleatorio correspondiente, y la matriz de covarianza de los efectos aleatorios es diagonal, es decir, efectos aleatorios no correlacionados.nlmefitsa

El,, y otros valores que devuelve esta función son el resultado de una simulación aleatoria (Montecarlo) diseñada para converger a las estimaciones de máxima verosimilitud de los parámetros.BETAPSI Debido a que los resultados son aleatorios, es aconsejable examinar la trama de simulación a los resultados para asegurarse de que la simulación ha convergido. También puede ser útil ejecutar la función varias veces, utilizando varios valores iniciales, o usar el parámetro para realizar varias simulaciones.'Replicates'

[BETA,PSI,STATS,B] = nlmefitsa(X,Y,GROUP,V,MODELFUN,BETA0,'Name',Value) acepta uno o más pares de nombre/valor de parámetro separados por comas. Especificar Name dentro de comillas simples.

Argumentos de entrada

Definitions:

En la siguiente lista de argumentos, se aplican las siguientes definiciones de variables:

  • — el número de observacionesn

  • — número de variables predictorash

  • — número de gruposm

  • — número de variables predictoras específicas del grupog

  • — el número de parámetrosp

  • — número de efectos fijosf

X

Una-por-matriz de observaciones sobre variables predictoras.nhnh

Y

Un vector de respuestas de un-por-1.n

GROUP

Una variable de agrupación que indica a cuál de los grupos pertenece cada observación. puede ser una variable categórica, un vector numérico, una matriz de caracteres con filas para los nombres de grupo, una matriz de cadenas o una matriz de celdas de vectores de caracteres.mGROUP

V

Una por matriz de variables predictoras específicas del grupo para cada uno de los grupos de los datos.mggm Estos son valores predictores que toman el mismo valor para todas las observaciones de un grupo. Las filas de se ordenan según.VGRP2IDX(GROUP) Utilice una matriz por celda para si cualquiera de los valores predictores específicos del grupo varía en tamaño entre grupos.mgV Especifique si no hay predictores de grupo.[]V

MODELFUN

Un identificador de una función que acepta valores predictores y parámetros de modelo, y devuelve valores ajustados. tiene el formulario con argumentos de entradaMODELFUNYFIT = MODELFUN(PHI,XFUN,VFUN)

  • — 1 por vector de los parámetros del modelo.PHIp

  • — Una-por-matriz de variables predictoras dondeXFUNlh

    • es 1 si es una sola fila delXFUNX

    • Esl Ni Si contiene las filas de un solo grupo de tamañoXFUNX Ni

    • es si contiene todas las filas de.lnXFUNX

  • — CualquieraVFUN

    • 1 por vector de predictores específicos de grupo para un solo grupo, correspondiente a una sola fila degV

    • Una-por-matriz, donde la-ésima fila de es (,:) Si la observación-TH está en grupo.ngkVFUNViki

    Si está vacío, las llamadas con sólo dos entradas.VnlmefitsaMODELFUN

Devuelve un vector-by-1 de valores ajustados.MODELFUNlYFIT Cuando una o contiene una sola fila, esa fila corresponde a todas las filas de los otros dos argumentos de entrada.PHIVFUN Para mejorar el rendimiento, utilice el par nombre/valor de parámetro (descrito a continuación) si puede calcular para más de un vector de parámetros de modelo en una llamada.'Vectorization'MODELFUNYFIT

BETA0

Vector a-by-1 con estimaciones iniciales para los efectos fijos.ff De forma predeterminada, es igual al número de parámetros del modelo. también puede ser un-por-matriz, y la estimación se repite veces utilizando cada columna de como un conjunto de valores iniciales.f pBETA0fREPSREPSBETA0

Argumentos de par nombre-valor

De forma predeterminada, se ajusta a un modelo donde cada parámetro de modelo es la suma de un efecto fijo y aleatorio correspondiente.nlmefitsa Utilice los siguientes pares de nombre/valor de parámetro para ajustar un modelo con un número diferente de o dependencia de efectos fijos o aleatorios. Utilice como máximo un nombre de parámetro con un prefijo y un nombre de parámetro con un prefijo.'FE''RE' Tenga en cuenta que algunas opciones cambian la forma de llamadas, como se describe a continuación.nlmefitsaMODELFUN

'FEParamsSelect'

Un vector que especifica qué elementos del vector de parámetro de modelo incluyen un efecto fijo, como un vector numérico con elementos en 1:, o como un vector 1 por lógica.PHIpp El modelo incluirá efectos fijos, donde es el número especificado de elementos.ff

'FEConstDesign'

Matriz A-by-Design, donde * son los componentes fijos de los elementos de.pfADESIGNADESIGNBETApPHI

'FEGroupDesign'

A-by--by-array especificando una matriz de diseño de efectos diferente por fijo para cada uno de los grupos.pfmpfm

'REParamsSelect'

Un vector que especifica qué elementos del vector de parámetro de modelo incluyen un efecto aleatorio, como un vector numérico con elementos en 1:, o como un vector 1 por lógica.PHIpp El modelo incluirá efectos aleatorios, donde es el número especificado de elementos.rr

'REConstDesign'

Matriz A-by-Design, donde * son los componentes aleatorios de los elementos de.prBDESIGNBDESIGNBpPHI Esta matriz debe constar de 0s y 1s, con un máximo de 1 por fila.

El modelo predeterminado es equivalente a la configuración de ambos y a, o para establecer ambos y a 1:.FEConstDesignREConstDesigneye(p)FEParamsSelectREParamsSelectp

Los pares adicionales de nombre/valor de parámetro opcionales controlan el algoritmo iterativo utilizado para maximizar la probabilidad:

'CovPattern'

Especifica una matriz numérica o por lógica que define el patrón de la matriz de covarianza de efectos aleatorios. calcula estimaciones para las desviaciones a lo largo de la diagonal, así como las covarianzas que corresponden a ceros en la diagonal de. restringe las covarianzas restantes, es decir, las que corresponden a ceros fuera de la diagonal, para que sean cero. debe ser una permutación de columna de fila de una matriz diagonal de bloque y agrega elementos que no sean cero según sea necesario para producir este tipo de patrón.rrPATPSInlmefitsaPSIPATnlmefitsaPATPATnlmefitsaPAT El valor predeterminado es, correspondiente a los efectos aleatorios no correlacionados.PATeye(r)

Como alternativa, especifique como un 1 por vector que contiene valores en 1:.PATrr En este caso, los elementos con valores iguales definen grupos de efectos aleatorios, calcula las covarianzas solo dentro de los grupos y restringe las covarianzas entre los grupos para que sean cero.PATnlmefitsa

'Cov0'

Valor inicial de la matriz de covarianza.PSI Debe ser una matriz definida por positivo.rr Si está vacío, el valor predeterminado depende de los valores de.BETA0

'ComputeStdErrors'

para calcular los errores estándar para las estimaciones de coeficiente y almacenarlos en la estructura de salida, o (predeterminado) para omitir este cálculo.trueSTATSfalse

'ErrorModel'

Un vector de caracteres o un escalar de cadena que especifica la forma del término de error. El valor predeterminado es.'constant' Cada modelo define el error utilizando una variable normal estándar (gaussiana), el valor de la función y uno o dos parámetros y.efab Las opciones son

  • 'constant' y = f + a*e

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

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

  • 'exponential' y = f*exp(a*e), o equivalentemente log(y) = log(f) + a*e

Si se da este parámetro, el campo de salida tiene el valorSTATS.errorparam

  • para ya'constant''exponential'

  • Parab'proportional'

  • [] paraab'combined'

'ErrorParameters'

Un vector escalar o de dos elementos que especifica los valores iniciales para los parámetros del modelo de error. Esto especifica los valores, o [] dependiendo del parámetro.ababErrorModel

'LogLikMethod'

Especifica el método para aproximar el logverosimilitud. Las opciones son:

  • — Muestreo de importancia'is'

  • — Cuadratura gaussiana'gq'

  • — Linearización'lin'

  • — Omitir la aproximación de logverosimilitud (por defecto)'none'

'NBurnIn'

Número de iteraciones de grabación iniciales durante las cuales no se recalculan las estimaciones de parámetros. El valor predeterminado es 5.

'NChains'

Número de "cadenas" simuladas.c El valor predeterminado es 1. Establecer > 1 hace que se calculen vectores de coeficiente simulados para cada grupo durante cada iteración.cc El valor predeterminado depende de los datos y se elige para proporcionar alrededor de 100 grupos en todas las cadenas.

'NIterations'

Número de iteraciones. Puede tratarse de un vector escalar o de tres elementos. Controla cuántas iteraciones se realizan para cada una de las tres fases del algoritmo:

  1. de recocido simulado

  2. tamaño de paso completo

  3. tamaño de paso reducido

El valor predeterminado es.[150 150 100] Un escalar se distribuye A través de las tres fases en las mismas proporciones que el valor predeterminado.

'NMCMCIterations'

Número de iteraciones de la cadena Markov Monte Carlo (MCMC). Puede tratarse de un vector escalar o de tres elementos. Controla cuántos de los tres tipos diferentes de actualizaciones MCMC se realizan durante cada fase de la iteración principal:

  1. actualización multivariada completa

  2. única actualización de coordenadas

  3. actualización de varias coordenadas

El valor predeterminado es.[2 2 2] Un valor escalar se trata como un vector de tres elementos con todos los valores iguales al escalar.

'OptimFun'

O bien, especificando la función de optimización que se utilizará durante el proceso de estimación.'fminsearch''fminunc' El valor predeterminado es.'fminsearch' Uso de requiere.'fminunc'Optimization Toolbox™

'Options'

Una estructura creada por una llamada a. utiliza los siguientes parámetros:statsetnlmefitsastatset

  • — Diferencia relativa utilizada en el cálculo del gradiente de diferencia finita.'DerivStep' Puede ser un escalar, o un vector cuya longitud es el número de parámetros del modelo.p El valor predeterminado es.eps^(1/3)

  • — Nivel de visualización durante la estimación.Display

    • (predeterminado): no muestra información'off'

    • : Muestra información después de la iteración final del algoritmo de estimación'final'

    • : Muestra información en cada iteración'iter'

  • FunValCheck

    • (sdefault): comprueba si hay valores no válidos (como o) de'on'NaNInfMODELFUN

    • — Omita esta comprobación'off'

  • — Identificador de función especificado mediante, una matriz de celdas con identificadores de función o una matriz vacía. llama a todas las funciones de salida después de cada iteración.OutputFcn@nlmefitsa Consulte (la función de salida predeterminada) para obtener un ejemplo de una función de salida.nlmefitoutputfcn.mnlmefitsa

'ParamTransform'

Un vector de-Values que especifica una función de transformación para cada uno de los parámetros:pf()p Cada elemento del vector debe ser uno de los siguientes códigos enteros especificando la transformación para el valor correspondiente de:

XB = ADESIGN*BETA + BDESIGN*B  PHI = f(XB) 
PHI

  • 0: (por defecto para todos los parámetros)PHI = XB

  • 1:log(PHI) = XB

  • 2:probit(PHI) = XB

  • 3:logit(PHI) = XB

'Replicates'

Número de estimaciones que se realizarán a partir de los valores iniciales del vector.REPSBETA0 Si es una matriz, debe coincidir con el número de columnas.BETA0REPSBETA0 El valor predeterminado es el número de columnas en.BETA0

'Vectorization'

Determina los posibles tamaños de los argumentos de entrada y, a.PHIXFUNVFUNMODELFUN Los valores posibles son:

  • : es una función (como un solucionador de ODE) que sólo puede calcular para un único conjunto de parámetros de modelo a la vez, es decir, debe ser un vector de fila única en cada llamada. llamadas en un bucle si es necesario utilizando un solo vector y con filas contenedora para una sola observación o grupo a la vez. puede ser una sola fila que se aplica a todas las filas de, o una matriz con filas correspondientes a filas en.'SinglePhi'MODELFUNYFITPHInlmefitsaMODELFUNPHIXFUNVFUNXFUNXFUN

  • — sólo puede aceptar entradas correspondientes a un solo grupo en los datos, es decir, debe contener filas de un solo grupo en cada llamada.'SingleGroup'MODELFUNXFUNX Dependiendo del modelo, es una sola fila que se aplica a todo el grupo, o una matriz con una fila para cada observación. es una sola fila.PHIVFUN

  • : puede aceptar entradas para varios vectores de parámetros y varios grupos en los datos.'Full'MODELFUN O puede ser una sola fila que se aplica a todas las filas de, o una matriz con filas correspondientes a filas en.PHIVFUNXFUNXFUN El uso de esta opción puede mejorar el rendimiento reduciendo el número de llamadas a, pero puede requerir para realizar la expansión de singleton en o.MODELFUNMODELFUNPHIV

El valor predeterminado es.'Vectorization''SinglePhi' En todos los casos, si está vacío, las llamadas con sólo dos entradas.VnlmefitsaMODELFUN

Argumentos de salida

BETA

Las estimaciones de los efectos fijos

PSI

Matriz de covarianza de una estimación para los efectos aleatorios.rr De forma predeterminada, es igual al número de parámetros del modelo.rp

STATS

Una estructura con los siguientes campos:

  • — La logverosimilitud maximizada para el modelo ajustado; vacía si el parámetro tiene su valor predeterminado deloglLogLikMethod'none'

  • — La raíz cuadrada de la varianza de error estimada (calculada en la escala de registro para el modelo de error)rmseexponential

  • — Los parámetros estimados del modelo de desviación de errorerrorparam

  • — El criterio de información de Akaike (vacío si está vacío), calculado como = – 2 * + 2 *, dondeaicloglaicloglnumParam

    • es la máxima logverosimilitud.logl

    • es el número de parámetros de ajuste, incluido el grado de libertad para la matriz de covarianza de los efectos aleatorios, el número de efectos fijos y el número de parámetros del modelo de error.numParam

  • — El criterio de información bayesiana (vacío si está vacío), calculado como =-2 * + log () *bicloglbicloglMnumParam

    • es el número de grupos.M

    • y se definen como in.loglnumParamaic

    Tenga en cuenta que cierta literatura sugiere que el cómputo de debe ser, =-2 * + log () *, donde está el número de observaciones.bicbicloglNnumParamN Para ajustar el valor de la salida se puede redefinir de la siguiente manera: =-(()) + ()bicbicbicnumeluniquegroupnumelY

  • — Los errores estándar para BETA (vacíos si el parámetro tiene su valor predeterminado de false)sebetaComputeStdErrors

  • — La covarianza estimada de las estimaciones de parámetros (vacía si es falsa)covbComputeStdErrors

  • — Los grados de error de libertaddfe

  • — Los residuos de población, ¿dónde están los valores previstos de la poblaciónpres(y-y_population)y_population

  • — Los residuales de población, donde se encuentran los valores previstos individualesires(y-y_population)y_population

  • — Los residuos ponderados de la poblaciónpwres

  • — Los residuales ponderados condicionalescwres

  • — Los residuos ponderados individualesiwres

Ejemplos

contraer todo

Cargue los datos de ejemplo.

load indomethacin 

Ajustar un modelo a los datos sobre las concentraciones de la droga indometacina en el torrente sanguíneo de seis sujetos durante ocho horas.

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  

Trace los datos junto con un ajuste de población general

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

Trazar curvas individuales basadas en estimaciones de efecto aleatorio.

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 

Algoritmos

Con el fin de estimar los parámetros de un modelo de efectos mixtos no lineales, nos gustaría elegir los valores de parámetro que maximizan una función de probabilidad. Estos valores se denominan estimaciones de máxima verosimilitud. La función de probabilidad puede escribirse en la forma

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

Dónde

  • son los datos de respuestay

  • β es el vector de los coeficientes de población

  • Σ2 es la desviación residual

  • ∑ es la matriz de covarianza para los efectos aleatorios

  • es el conjunto de efectos aleatorios no observadosb

Cada () función en el lado derecho es una función de probabilidad normal (gaussiana) que puede depender de covariables.p

Dado que la integral no tiene una forma cerrada, es difícil encontrar parámetros que lo maximicen. DeLyon, Lavielle y Moulines propusieron encontrar las estimaciones de máxima verosimilitud utilizando un algoritmo de maximización de expectativas (EM) en el que el paso E se sustituye por un procedimiento estocástico.[1] Llamaron a su algoritmo SAEM, para la aproximación estocástica EM. Demostraron que este algoritmo tiene propiedades teóricas deseables, incluyendo la convergencia en condiciones prácticas y la convergencia a un máximo local de la función de probabilidad. Su propuesta involucra tres pasos:

  1. Simulación: Genere valores simulados de los efectos aleatorios de la densidad posterior (| Σ) dadas las estimaciones de parámetros actuales.bpb

  2. Aproximación estocástica: Actualice el valor esperado de la función logverosimilitud tomando su valor del paso anterior y moviendo la parte hacia el valor medio de la logverosimilitud calculada a partir de los efectos aleatorios simulados.

  3. Paso de maximización: Elija nuevas estimaciones de parámetros para maximizar la función logverosimilitud dados los valores simulados de los efectos aleatorios.

Referencias

[1] Delyon, B., M. Lavielle, and E. Moulines, Convergence of a stochastic approximation version of the EM algorithm, Annals of Statistics, 27, 94-128, 1999.

[2] Mentré, F., and M. Lavielle, Stochastic EM algorithms in population PKPD analyses, American Conference on Pharmacometrics, 2008.

Introducido en R2010a