# Linear regression with a given covariance structure in R

I want to fit a linear model in R with a given covariance structure: $$y=X\beta+\epsilon$$ where the covariance matrix of $$\epsilon$$ is block diagonal by a grouping factor. Suppose there are $$B$$ blocks. Within block $$b$$, the covariance is of the following form: $$\sigma_1^211^T+\sigma_2^2G_b+\sigma_3^2I$$ where $$1$$ is the vector of all ones, $$I$$ is the identity matrix and $$G_b$$ is generated by a Gaussian process kernel: $$G_{b}(i,j)=k(x_{bi},x_{bj},\theta)$$. Here $$x_{bi}$$ and $$x_{bj}$$ are sample $$i,j$$ in block $$b$$. $$\theta$$ is an unknown parameter of the kernel $$k$$. $$\sigma_1,\sigma_2,\sigma_3,\theta$$ are all unknown and are the same across all blocks. An example of $$k$$ is $$k(x_{bi},x_{bj},\theta)=\exp(-\frac{||x_{bi}-x_{bj}||^2}{\theta^2})$$

Is there any existing packages in R to fit the above model and implement inference?

• It is easy to fit this model in SAS by PROC MIXED. Nov 8, 2018 at 3:49

You could fit this type of model using the lme() function from the nlme package; a generic example is:
lme(y ~ x1 + x2 + ..., data = <your_data>, random = ~ 1 | block_id,

The only difference from your postulated covariance matrix will be that in the above code you have $$\sigma^2 = \sigma_2^2 + \sigma_3^2$$.