Hi R users,

I need to write a function that considers the row and column 
autoregressive residual correlation structures. Usually 
R=sigma^2_e*Identity matrix of order n, where n is the number of 
observations and sigma^2_e is the usual residual error. I have to 
consider an AR1 * AR1 where * stands for the Kronecker product. The 
formula  is R=sigma^2_e times Vc(phi_c) * Vr(phi_r) where
* is Kronecker product as explained above, Vc(phi_c) is the AR1 
correlated structure (matrix) on the columns (y coordinate) of order n*n 
and Vc(phi_r) is the AR1 correlated structure (matrix) on the rows (x 
coordinate) of order n*n with phi_r and phi_c denoting  the spatial 
correlation value in the AR1 structure.

For example. If I have a field experimental design called des1 with x 
and y coordinates given, number of blocks and the genotypes applied to 
each of these blocks:

des1
   x y block genotypes
1 1 1     1         4
2 2 1     1         3
3 1 2     1         1
4 2 2     1         2
5 3 1     2         4
6 4 1     2         2
7 3 2     2         3
8 4 2     2         1

If I  assume that sigma^2_e =1, phi_c=phi_r = 0.1, then the Residual structure should be

[,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
[1,]  1.02030 -0.10203 -0.10203  0.01020  0.00000  0.00000  0.00000  0.00000
[2,] -0.10203  1.03051  0.01020 -0.10305 -0.10203  0.00000  0.01020  0.00000
[3,] -0.10203  0.01020  1.02030 -0.10203  0.00000  0.00000  0.00000  0.00000
[4,]  0.01020 -0.10305 -0.10203  1.03051  0.01020  0.00000 -0.10203  0.00000
[5,]  0.00000 -0.10203  0.00000  0.01020  1.03051 -0.10203 -0.10305  0.01020
[6,]  0.00000  0.00000  0.00000  0.00000 -0.10203  1.02030  0.01020 -0.10203
[7,]  0.00000  0.01020  0.00000 -0.10203 -0.10305  0.01020  1.03051 -0.10203
[8,]  0.00000  0.00000  0.00000  0.00000  0.01020 -0.10203 -0.10203  1.02030


I wrote an R function which works well and takes a second or less to give me the above matrix. However, it takes too long time (about 30 minutes) when I have very large dataset (e.g. with n=5000 observations).
Here is the code I wrote.

# An R function to calculate R and R_inverse from spatial data
rspat<-function(rhox,rhoy,s2e=1)
{
require matlab
   R<-s2e*eye(N) # library(matlab required)
   for(i in 1:N) {
     for (j in i:N){
       y1<-y[i]
       y2<-y[j]
       x1<-x[i]
       x2<-x[j]
       R[i,j]<-s2e*(rhox^abs(x2-x1))*(rhoy^abs(y2-y1)) # Core AR(1)*AR(1)
       R[j,i]<-R[i,j]
     }
   }
   IR<-ginv(R) # use the Penrose Generalized inverse
   IR
}

Is there an existing/in built R code that does the same in a more 
efficient and faster than my code?

Thanks

	[[alternative HTML version deleted]]

