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.

Modelos de efectos mixtos con nlmefit y nlmefitsa

Cargue los datos de ejemplo.

load indomethacin 

Los datos en registros concentraciones de la droga indometacina en el torrente sanguíneo de seis sujetos durante ocho horas.indomethacin.mat

Trazar el diagrama de dispersión de la indometacina en el torrente sanguíneo agrupado por tema.

gscatter(time,concentration,subject) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') title('{\bf Indomethacin Elimination}') hold on 

página se describe un modelo útil para este tipo de datos.Especificación de modelos de efectos mixtos

Construya el modelo a través de una función anónima.

model = @(phi,t)(phi(1)*exp(-exp(phi(2))*t) + ...                  phi(3)*exp(-exp(phi(4))*t)); 

Utilice la función para ajustar el modelo a todos los datos, ignorando los efectos específicos del sujeto.nlinfit

phi0 = [1 2 1 1]; [phi,res] = nlinfit(time,concentration,model,phi0); 

Calcule el error medio cuadrado.

numObs = length(time); numParams = 4; df = numObs-numParams; mse = (res'*res)/df 
 mse =      0.0304  

Súper imponer el modelo en el gráfico de dispersión de datos.

tplot = 0:0.01:8; plot(tplot,model(phi,tplot),'k','LineWidth',2) hold off 

Dibuje el diagrama de caja de los residuos por tema.

colors = 'rygcbm'; h = boxplot(res,subject,'colors',colors,'symbol','o'); set(h(~isnan(h)),'LineWidth',2) hold on boxplot(res,subject,'colors','k','symbol','ko') grid on xlabel('Subject') ylabel('Residual') hold off 

El diagrama de caja de los residuos por asunto muestra que las casillas están sobre todo por encima o por debajo de cero, lo que indica que el modelo no ha tenido en cuenta los efectos específicos del sujeto.

Para tener en cuenta los efectos específicos del sujeto, ajuste el modelo por separado a los datos para cada tema.

phi0 = [1 2 1 1]; PHI = zeros(4,6); RES = zeros(11,6); for I = 1:6     tI = time(subject == I);     cI = concentration(subject == I);     [PHI(:,I),RES(:,I)] = nlinfit(tI,cI,model,phi0); end PHI 
 PHI =      2.0293    2.8277    5.4683    2.1981    3.5661    3.0023     0.5794    0.8013    1.7498    0.2423    1.0408    1.0882     0.1915    0.4989    1.6757    0.2545    0.2915    0.9685    -1.7878   -1.6354   -0.4122   -1.6026   -1.5069   -0.8731  

Calcule el error medio cuadrado.

numParams = 24; df = numObs-numParams; mse = (RES(:)'*RES(:))/df 
 mse =      0.0057  

Trace el gráfico de dispersión de los datos y superponga el modelo para cada asignatura.

gscatter(time,concentration,subject) xlabel('Time (hours)') ylabel('Concentration (mcg/ml)') title('{\bf Indomethacin Elimination}') hold on for I = 1:6     plot(tplot,model(PHI(:,I),tplot),'Color',colors(I)) end axis([0 8 0 3.5]) hold off 

da estimaciones de los cuatro parámetros del modelo para cada uno de los seis sujetos.PHI Las estimaciones varían considerablemente, pero se toman como un modelo de 24 parámetros de los datos, el error cuadrático medio de 0,0057 es una reducción significativa de 0,0304 en el modelo de cuatro parámetros original.

Dibuje el diagrama de caja de los residuos por tema.

h = boxplot(RES,'colors',colors,'symbol','o'); set(h(~isnan(h)),'LineWidth',2) hold on boxplot(RES,'colors','k','symbol','ko') grid on xlabel('Subject') ylabel('Residual') hold off 

Ahora la gráfica de caja muestra que el modelo más grande representa la mayoría de los efectos específicos del sujeto. La propagación de los residuos (la escala vertical del diagrama de caja) es mucho menor que en la gráfica de cuadro anterior, y las casillas ahora se centran principalmente en cero.

Mientras que el modelo de 24 parámetros cuenta con éxito las variaciones debido a los temas específicos en el estudio, no considera los temas como representantes de una población más grande. La distribución de muestreo de la que se dibujan los sujetos es probablemente más interesante que la muestra en sí. El propósito de los modelos de efectos mixtos es tener en cuenta las variaciones específicas del sujeto de manera más amplia, ya que los efectos aleatorios varían en función de la población.

Utilice la función para ajustar un modelo de efectos mixtos a los datos.nlmefit También se puede utilizar en lugar de.nlmefitsanlmefit

La siguiente función anónima,,, adapta el modelo de cuatro parámetros utilizado por la sintaxis de llamada al permitir que se separen los distintos elementos para cada individuo.nlme_modelnlinfitnlmefit De forma predeterminada, asigna efectos aleatorios a todos los parámetros del modelo.nlmefit También de forma predeterminada, asume una matriz de covarianza diagonal (sin covarianza entre los efectos aleatorios) para evitar la sobreparametrización y los problemas de convergencia relacionados.nlmefit

nlme_model = @(PHI,t)(PHI(:,1).*exp(-exp(PHI(:,2)).*t) + ...                       PHI(:,3).*exp(-exp(PHI(:,4)).*t)); phi0 = [1 2 1 1]; [phi,PSI,stats] = nlmefit(time,concentration,subject, ...                           [],nlme_model,phi0) 
 phi =      2.8277     0.7729     0.4606    -1.3459   PSI =      0.3264         0         0         0          0    0.0250         0         0          0         0    0.0124         0          0         0         0    0.0000   stats =     struct with fields:             dfe: 57           logl: 54.5882            mse: 0.0066           rmse: 0.0787     errorparam: 0.0815            aic: -91.1765            bic: -93.0506           covb: [4x4 double]         sebeta: [0.2558 0.1066 0.1092 0.2244]           ires: [66x1 double]           pres: [66x1 double]          iwres: [66x1 double]          pwres: [66x1 double]          cwres: [66x1 double]  

El error medio cuadrado de 0,0066 es comparable al 0,0057 del modelo de 24 parámetros sin efectos aleatorios, y significativamente mejor que el 0,0304 del modelo de cuatro parámetros sin efectos aleatorios.

La matriz de covarianza estimada muestra que la varianza del cuarto efecto aleatorio es esencialmente cero, lo que sugiere que se puede eliminar para simplificar el modelo.PSI Para ello, utilice el par nombre-valor para especificar los índices de los parámetros que se modelan con efectos aleatorios en.'REParamsSelect'nlmefit

[phi,PSI,stats] = nlmefit(time,concentration,subject, ...                           [],nlme_model,phi0, ...                           'REParamsSelect',[1 2 3]) 
 phi =      2.8277     0.7728     0.4605    -1.3460   PSI =      0.3270         0         0          0    0.0250         0          0         0    0.0124   stats =     struct with fields:             dfe: 58           logl: 54.5875            mse: 0.0066           rmse: 0.0780     errorparam: 0.0815            aic: -93.1750            bic: -94.8410           covb: [4x4 double]         sebeta: [0.2560 0.1066 0.1092 0.2244]           ires: [66x1 double]           pres: [66x1 double]          iwres: [66x1 double]          pwres: [66x1 double]          cwres: [66x1 double]  

La log-verosimilitud es casi idéntica a lo que era con efectos aleatorios para todos los parámetros, el criterio de información de Akaike se reduce de-91,1765 a-93,1750, y el criterio de información bayesiana se reduce de-93,0506 a-94,8410.loglaicbic Estas medidas apoyan la decisión de eliminar el cuarto efecto aleatorio.

La reposición del modelo simplificado con una matriz de covarianza completa permite la identificación de correlaciones entre los efectos aleatorios. Para ello, utilice el parámetro para especificar el patrón de elementos distintos de cero en la matriz de covarianza.CovPattern

[phi,PSI,stats] = nlmefit(time,concentration,subject, ...                           [],nlme_model,phi0, ...                           'REParamsSelect',[1 2 3], ...                           'CovPattern',ones(3)) 
 phi =      2.8148     0.8293     0.5613    -1.1407   PSI =      0.4767    0.1152    0.0499     0.1152    0.0321    0.0032     0.0499    0.0032    0.0236   stats =     struct with fields:             dfe: 55           logl: 58.4731            mse: 0.0061           rmse: 0.0782     errorparam: 0.0781            aic: -94.9463            bic: -97.2369           covb: [4x4 double]         sebeta: [0.3028 0.1103 0.1179 0.1662]           ires: [66x1 double]           pres: [66x1 double]          iwres: [66x1 double]          pwres: [66x1 double]          cwres: [66x1 double]  

La matriz de covarianza estimada muestra que los efectos aleatorios en los dos primeros parámetros tienen una correlación relativamente fuerte, y ambos tienen una correlación relativamente débil con el último efecto aleatorio.PSI Esta estructura en la matriz de covarianza es más evidente si se convierte a una matriz de correlación mediante.PSIcorrcov

RHO = corrcov(PSI) clf; imagesc(RHO) set(gca,'XTick',[1 2 3],'YTick',[1 2 3]) title('{\bf Random Effect Correlation}') h = colorbar; set(get(h,'YLabel'),'String','Correlation'); 
 RHO =      1.0000    0.9316    0.4707     0.9316    1.0000    0.1178     0.4707    0.1178    1.0000  

Incorpore esta estructura en el modelo cambiando la especificación del patrón de covarianza a Block-diagonal.

P = [1 1 0;1 1 0;0 0 1] % Covariance pattern [phi,PSI,stats,b] = nlmefit(time,concentration,subject, ...                             [],nlme_model,phi0, ...                             'REParamsSelect',[1 2 3], ...                             'CovPattern',P) 
 P =       1     1     0      1     1     0      0     0     1   phi =      2.7830     0.8981     0.6581    -1.0000   PSI =      0.5180    0.1069         0     0.1069    0.0221         0          0         0    0.0454   stats =     struct with fields:             dfe: 57           logl: 58.0804            mse: 0.0061           rmse: 0.0768     errorparam: 0.0782            aic: -98.1608            bic: -100.0350           covb: [4x4 double]         sebeta: [0.3171 0.1073 0.1384 0.1453]           ires: [66x1 double]           pres: [66x1 double]          iwres: [66x1 double]          pwres: [66x1 double]          cwres: [66x1 double]   b =     -0.8507   -0.1563    1.0427   -0.7559    0.5652    0.1550    -0.1756   -0.0323    0.2152   -0.1560    0.1167    0.0320    -0.2756    0.0519    0.2620    0.1064   -0.2835    0.1389  

La estructura de covarianza en diagonal de bloque se reduce de-94,9462 a-98,1608 y de-97,2368 a-100,0350 sin afectar significativamente al log-verosimilitud.aicbic Estas medidas soportan la estructura de covarianza utilizada en el modelo final. La salida da predicciones de los tres efectos aleatorios para cada uno de los seis sujetos.b Estos se combinan con las estimaciones de los efectos fijos en para producir el modelo de efectos mixtos.phi

Trazar el modelo de efectos mixtos para cada uno de los seis sujetos. Para la comparación, también se muestra el modelo sin efectos aleatorios.

PHI = repmat(phi,1,6) + ...                 % Fixed effects       [b(1,:);b(2,:);b(3,:);zeros(1,6)];    % Random effects RES = zeros(11,6); % Residuals colors = 'rygcbm'; for I = 1:6     fitted_model = @(t)(PHI(1,I)*exp(-exp(PHI(2,I))*t) + ...                         PHI(3,I)*exp(-exp(PHI(4,I))*t));     tI = time(subject == I);     cI = concentration(subject == I);     RES(:,I) = cI - fitted_model(tI);      subplot(2,3,I)     scatter(tI,cI,20,colors(I),'filled')     hold on     plot(tplot,fitted_model(tplot),'Color',colors(I))     plot(tplot,model(phi,tplot),'k')     axis([0 8 0 3.5])     xlabel('Time (hours)')     ylabel('Concentration (mcg/ml)')     legend(num2str(I),'Subject','Fixed') end 

Si se omiten los valores atípicos obvios en los datos (visibles en los trazados de cuadro anteriores), una gráfica de probabilidad normal de los residuos muestra un acuerdo razonable con los supuestos del modelo sobre los errores.

clf; normplot(RES(:))