# Generalized linear regression with custom variance-covariance matrix in R

I want to compute the estimate of $$\beta$$ for a linear model $$Y = X\beta + \varepsilon$$ with $$\varepsilon \sim N_d(0, \sigma^2V),$$ where $$V$$ is a $$d\times d$$ definitive posive, symmetric matrix.

It is straightforward to generalize the theory of ordinary linear models (where $$\varepsilon \sim N_d(0, \sigma^2I_d))$$, to this general case. In fact, it is sufficient to consider the Cholesky decomposition of the matrix $$V$$ and then transform the variables (details, for instance, here).

Now that the theory is clear, how can I specify and solve this problem in R?

• You can use the gls function in the nlme library cran.r-project.org/web/packages/nlme . Technically this is called generalized least squares (GLS), the errors are unequal or correlated but still normal. In generalized linear models (GLM), the errors follow some distribution other than the normal distribution. – StatGrrl Feb 10 at 12:57
• Thank you for the package suggestion and the clarification of the terminology – Nisba Feb 10 at 13:16
• If V is a diagonal matrix then the errors have unequal variances but are independent - this is weighted least squares (WLS).The gls function allows for adding a correlation structure and/or weights. But if all you need is WLS, the lm function allows for this. – StatGrrl Feb 10 at 13:18
• I need to consider a correlation structure – Nisba Feb 10 at 13:21
• Is V known or must it be estimated? – StatGrrl Feb 10 at 13:31