I'm working on a project where I try to determine the influence of consulting expenditures companies made the year before on several economic figures. My dataset consists of roughly 2'000 different companies together with data for 2-7 years, thereby all continuous variables are log-transformed and vary roughly on the same scale.
To account for the occuring correlations I use a linear mixed-effect model (R lme) with the intercept as the only random effect. Under this setting, strong heteroscedasticity occurs and since I'm unable to determine any variance covariate, I wish to use
weights=varIdent(form=~1|Subject) to assign each company its own variance. However, using the
weights command in this way leads to a non-converging model.
To get convergence, I tried two things:
- I tried to estimate the residual variance by the original variance of the dependent variable and added
weights=varPower(form=~original_variance). While this leads to nicer looking residuals, the estimation produced obviously wrong results.
- I estimated the model first without the
weightsargument and used the resulting residual variance in a second step again as the argument to
weights=varPower(form=~residual_variance). This leads to better-looking residuals and at the same time the result seem to make a lot more sense.
Now my question: is approach 2 valid or do I bias the model in any way or produce wrong within-group variance-estimations and therefore wrong confidence intervals?