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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
