¿Cómo obtener/considerar la variabilidad para distintas variables categóricas en un modelo lineal de efectos mixtos?

Buenas, me encuentro realizando un modelo lineal de efectos mixtos de un conjunto de datos en donde tengo como efecto aleatorio un edificio, el cual se tienen distintos datos para un mismo edificio, cómo la altura, periodo, amortiguamiento, etc. Uno de ellos es el material del edificio considerada como variable categórica. Realicé el siguiente código:

lme = fitlme(Tabla_datos, 'logxi ~ logH * Material + (1|key) + (1|key:Material)', ...

'FitMethod', 'REML', ...

'CovariancePattern', {'Full', 'Diagonal'})

Donde key corresponde al nombre del edificio, este código me entrega los parámetros de covarianzas para el edificio, pero lo que yo necesito saber es cómo obtener la variabilidad para cada material en donde se trate a cada material como un estrato independiente con su propia variabilidad

Answer 1

[Umar]

Hi @Sofía Mendoza,

Sorry, I had to use Google translate to provide response to your comments. Instead of modeling the interaction between the building key and material as a single random effect, you can create separate random effects for each material. This would allow you to capture the variability in responses attributable to different materials independently. Here is how you might modify your model:

   lme = fitlme(data_table, 'logxi ~ logH + Material + (1|Material:key)', ...
                'FitMethod', 'REML', ...
                'CovariancePattern', {'Full', 'Diagonal'});

In this modified model, 1|Material:key allows you to estimate a random intercept for each combination of building material and building key, thus treating each material's influence independently while still accounting for the nested structure of your data. After fitting your model, you can extract the variance components associated with each random effect. This can be done using the `anova` function or by accessing the VarianceComponents property of your fitted LMM object:

   varComp = lme.VarianceComponents;

This will give you insight into how much variance is explained by the building and material interactions. To further validate your approach, consider comparing models with different structures using Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC). This will help you assess whether separating out materials leads to a significantly better fit compared to your original model. After obtaining your variance estimates, visualizing them can provide additional insights. You might create plots showing the distribution of logxi across different materials or boxplots comparing variability in responses across material types.

Make sure that your data meets the assumptions required for LMMs, including normality and homoscedasticity of residuals and consider using cross-validation techniques to assess the robustness of your findings across different subsets of your data.

Hope this helps.

[Sofía Mendoza]

Thank u so much for anwer!!

Disculpa te escribiré en español ya que no sé como expresarme en inglés (Chat GPT traduce de mejor manera que google). Estos son los resultados que estoy replicando:

Ya obtuve los valores para los efectos fijos, tanto coeficiente, error estándar etc. Sin embargo, aún no logro llegar a los resultados de variabilidad por material, ya que tengo entendido que MATLAB asume que en los efectos aleatorios utilizan un mismo sigma para todos, es por eso que quiero calcular esta variabilidad para las distintas categorías (en este caso: Materiales)

disculpa las molestias, espero puedas responder, saludos!

[Umar]

Hi @@Sofía Mendoza,

To calculate the variability by material, you can utilize the coefficients and standard errors provided in your regression output. The idea is to treat each material's effect as a random variable with its own variance, allowing us to estimate the variability across different materials. Here is how you can approach this problem using MATLAB. You will create a dataset representing your materials and building heights. In this case, I used synthetic data for it. Then, use MATLAB's fitlme function to fit a linear mixed-effects model.

https://www.mathworks.com/help/stats/fitlme.html

https://youtu.be/997ni1xcmKw?si=GknGYvSbFdbLvsyq

After fitting, extract the estimated variance components for each material category. Below is an illustrative example of how you might implement this in MATLAB:

% Generate synthetic data
numBuildings = 100; % Number of buildings
height = rand(numBuildings, 1) * 10; % Random heights between 0 
and 10   
meters
materials = randi([1, 5], numBuildings, 1); % Randomly assign 
materials (1-5)

% Create response variable based on fixed effects
fixedEffectIntercept = -1.84;
fixedEffectHeight = -0.38 * log(height);
materialEffects = [0; 0.52; 0.99; 0.32; -0.67]; % B2 to B6 
coefficients

% Generate response variable Y
Y = fixedEffectIntercept + fixedEffectHeight + 
materialEffects(materials) +    
randn(numBuildings, 1) * 0.5;

% Create a table for fitting
tbl = table(Y, height, materials);

% Fit mixed-effects model
lme = fitlme(tbl, 'Y ~ height + materials + (1|materials)');

% Display results
disp(lme);

% Extract variance components
varComponents = lme.Coefficients.Estimate;
varianceByMaterial = varComponents(2:end).^2; % Square of 
estimates gives 
variance

% Display final results
disp('Variability by Material:');
materialNames = {'Concrete', 'Steel', 'Mixed', 'Masonry', 
'Wood'};
disp(array2table(varianceByMaterial', 'VariableNames',   
materialNames(1:length(varianceByMaterial))));

Please see attached.

Once you have executed the code above and obtained the variance estimates for each material type, you can interpret these values as measures of how much variability exists within each category relative to the overall model. Make sure that your synthetic data reasonably reflects real-world conditions. The random noise added (randn(numBuildings, 1) * 0.5) should be adjusted based on expected error levels in your actual data.After fitting the model and obtaining variability estimates, consider conducting residual diagnostics or cross-validation to assess the robustness of your model.

Hope this helps.

[Sofía Mendoza]

Hola, agradezco mucho de tu ayuda, pero no logro llegar a los mismos resultador de la tabla. Quiero pensar que es debido a que en ningun momento consideraste logxi, pero aunque trate de incorporarlo no llego al resultado ya que el estudio es acerca de cómo afecta el material del edificio al amortiguamiento de este (logxi). Es por esta razón que al utilizar el siguiente código sí me dan los mismos resultados azules de la tabla (fixed effects).

lme = fitlme(Tabla_datos, 'logxi ~ logH * Material + (1|key) + (1|key:Material)', ...
             'FitMethod', 'REML', ...
             'CovariancePattern', {'Full', 'Diagonal'})

Pero no logro obtener los valores en amarillo de la tabla mostrada inicialmente :( ... Necesito obtener delta (δ) o sigma (σ) ya que posterior a eso puedo utilizar la siguiente ecuación:

Por otra parte quería aprovechar la oportunidad de preguntarte sobre esta parte del código que no comprendo muy bien cómo adaptarla a mis datos:

randn(numBuildings, 1) * 0.5;

Además, que función estaría cumpliendo la nueva función Y?

Disculpas las molestias, espero puedas ayudarme!!

[Torsten]

@Sofía Mendoza

Creo que podría obtener más respuestas a su pregunta si utilizara Google Translate español-inglés para sus preguntas antes de enviarlas.

[Sofía Mendoza]

Gracias!

[Umar]

Hi @Sofía Mendoza,

To address your concerns regarding the linear mixed-effects model and the specific results you are trying to achieve, let me break down the components of your model and the associated calculations step by step. The code snippet provided by you is to fit a linear mixed-effects model:

lme = fitlme(data_table, 'logxi ~ logH * Material + (1|key) + 
(1|key:Material)', ...'FitMethod', 'REML', ...
'CovariancePattern', {'Full', 'Diagonal'});

This model specifies that logxi (the response variable) is influenced by the interaction between logH and Material, while accounting for random effects associated with key and the interaction of key with Material. The key components of the model are

Fixed Effects: These are the coefficients that represent the average effect of the predictors on the response variable. In your output, you have estimates for the intercept, logH, and the different materials.

Random Effects: These account for variability in the data that is not explained by the fixed effects. The random effects structure you specified allows for different intercepts for each level of key and for each combination of key and Material.

Covariance Parameters: The covariance structure you chose (Full, Diagonal) affects how the random effects are modeled. A full covariance structure allows for correlations between random effects, while a diagonal structure assumes no correlation.

To obtain the values in yellow from the initial table, you need to ensure that your model is correctly specified and that the data is appropriately pre-processed. Here are some steps to consider:

Data Preparation: Ensure that your data_table is correctly formatted and that there are no missing values in the variables used in the model. Missing data can lead to unexpected results.

Model Specification: Double-check the interaction terms in your model. If you suspect that the interaction between logH and Material is not being captured correctly, consider simplifying the model or testing each term individually.

Calculating Delta (δ) or Sigma (σ): The equation you provided for calculating sigma involves the variances of the fixed effects and the residuals. You can extract these from the model output.

Now, let me address your query regarding clarification on random number generation

   randn(numBuildings, 1) * 0.5;

This generates a vector of random numbers from a standard normal distribution (mean = 0, standard deviation = 1) and scales them by 0.5. This is often used to simulate random effects or noise in your model. If you want to incorporate random effects into your model, you can use this to simulate variability among buildings.

The last query to address the role of function Y , you mentioned it likely represents the response variable in your model. In the context of mixed-effects models, Y would be the observed values of logxi that you are trying to predict based on the fixed and random effects specified in your model.

Hope this helps.