Exploring effects of different centering methods in hierarchical linear models on bias and standard errors through Monte-Carlo simulation
The hierarchical linear model described as follows:
The models were fitted using nlme package.
Two Level-1 predictors are generated randomly with means of 0.1 and 1.1 and variance of 1. Both Level-2 predictors have means of 1 and variance of 1.
The system is defined by ICC, number of groups, group size (all groups are assumed to have the same size), correlation between Level-1 predictors and list of gamma parameters.
Example. If all fixed effects coefficients are set to be 1, e.g. .
raw.data <- generate_data(group_size=10, num_groups=50, ICC=0.3,corr=0.7,gamma=rep(1,9))
High correlation between predictors may cause non-convergence. By default, generate_data function checks if the system is singular.
> set.seed(7)
> generate_data(group_size=10, num_groups=50, ICC=0.3,corr=0.7,gamma=rep(1,9))
[1] "Cgm model is not converging. Restarting simulation"
y x1 x2 g_id g1 g2
1 1.47476276 -0.2105077036 1.2572240503 1 -0.37154303 0.1137519
2 0.54211697 0.2740988978 1.8217880838 1 -0.37154303 0.1137519
3 1.23471334 -1.2305253710 0.4158803929 1 -0.37154303 0.1137519
...
Convergence check can be also done using conv_check function. If function returns 0, then convergence is not met.
> sample_data<-generate_data(group_size=10, num_groups=50, ICC=0.3,corr=0.7,gamma=rep(1,9))
> conv_check(sample_data)
[1] 1
cgm_cent(data) and cwc_cent(data) perform grand mean centering group mean centering of Level-1 predictors respectively.
> cwc_data<-cwc_cent(sample_data)
mlm_model(data) function is used to get a lme model object.
> raw_model<-raw_model(sample_data)
> cwc_model<-cwc_cent(sample_data)
To run a simulation, use simulation_run() function to generate raw data and get the parameters' estimates of raw scores, grand mean centered and group mean centered models. The results would be added to corresponding csv-files in /data
folder. Can be used together with replicate to get a desired number of observations. Returns 0 if fininshed sucessfully.
> replicate(n = 100, simulation_run(), simplify = FALSE)
To calculate bias, use get_coef_bias() function that takes two sets of parameters as an input (e.g. your raw coefficients and original gamma vector):
> coef.bias<-get_coef_bias(raw.coef,gamma)
Similar functions are avaliable for MSE and relative change:
> coef.mse<-get_coef_mse(raw.coef,gamma)
> cgm.relchange<-get_rel_change(raw.coef,cgm.coef)