[R] help in writing an R-function for Residual correlated structures

[ONKELINX, Thierry]

R works faster if you can avoid loops the loops. There is an example. Note that it required global variables (like your function). You better avoid that.

rspat <- function(rhox, rhoy, s2e = 1){
  require(matlab)
  R <- s2e * eye(N)
  i <- rep(seq_len(N), each = N)
  j <- rep(seq_len(N), N)
  j <- j[j > 1]
  R[cbind(i, j)] <- s2e*(rhox^abs(x[j] - x[i]))*(rhoy^abs(y[j] - y[i])) # Core AR(1)*AR(1)
  R[cbind(j, i)] <- R[cbind(i, j)]
  IR <- ginv(R) # use the Penrose Generalized inverse
  IR
}


ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be
www.inbo.be

To call in the statistician after the experiment is done may be no more than asking him to perform a post-mortem examination: he may be able to say what the experiment died of.
~ Sir Ronald Aylmer Fisher

The plural of anecdote is not data.
~ Roger Brinner

The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey


-----Oorspronkelijk bericht-----
Van: r-help-bounces op r-project.org [mailto:r-help-bounces op r-project.org] Namens Laz
Verzonden: donderdag 12 juni 2014 12:35
Aan: r-help op r-project.org
Onderwerp: [R] help in writing an R-function for Residual correlated structures

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 op 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.
* * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver weer en binden het INBO onder geen enkel beding, zolang dit bericht niet bevestigd is door een geldig ondertekend document.
The views expressed in this message and any annex are purely those of the writer and may not be regarded as stating an official position of INBO, as long as the message is not confirmed by a duly signed document.