2
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ It is easy to fit this model in SAS by PROC MIXED. $\endgroup$ – user158565 Nov 8 '18 at 3:49
2
$\begingroup$

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, 
    correlation = corGaus(form = ~ x1 | 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$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.