Laz noticed a slight discrepancy in my solution.  This can be fixed by
sorting the data first.  Here is a reproducible example using his original
method (for the non-inverted R matrix) and my version that does not use
(any obvious) loops.

Kevin Wright


des1 <-  data.frame(x=c(1,2,1,2,3,4,3,4), y=c(1,1,2,2,1,1,2,2),
                    block=c(1,1,1,1,2,2,2,2), gen=c(4,3,1,2,4,2,3,1))
des1 <- des1[order(des1$x),]
des1

rspat<-function(rhox, rhoy, s2e=1, N) {
  R <- s2e * diag(N)
  y <- des1$y
  x <- des1$x
  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]
    }
  }

  R
}
rspat(rhox = .6, rhoy = .6, s2e = 1, N = 8)

#
----------------------------------------------------------------------------

library(mvtnorm)
nr <- 4; nc <- 2
rho.r <- .6; rho.c <- 0.6
sigAR <- diag(nr)
sigAR <- rho.r^ abs(row(sigAR) - col(sigAR))
sigAC <- diag(nc)
sigAC <- rho.c ^ abs(row(sigAC) - col(sigAC))
sig <- 1 * kronecker(sigAR, sigAC)
round(sig,3)

# Compare
round(rspat(rhox = .6, rhoy = .6, s2e = 1, N = 8),3)
round(sig,3)




On Thu, Jun 12, 2014 at 11:52 PM, Laz <lmramba@ufl.edu> wrote:

>  I have switched them and the results is now very close to the original,
> with a few elements still interchanged as shown below:
>
> I expect the below matrix for rhox=rhoy=0.6
>
> r1          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
> [1,]  2.441406 -1.464844 -1.464844  0.878906  0.000000  0.000000  0.000000  0.000000
> [2,] -1.464844  3.320312  0.878906 -1.992187 -1.464844  0.000000  0.878906  0.000000
> [3,] -1.464844  0.878906  2.441406 -1.464844  0.000000  0.000000  0.000000  0.000000
> [4,]  0.878906 -1.992187 -1.464844  3.320312  0.878906  0.000000 -1.464844  0.000000
> [5,]  0.000000 -1.464844  0.000000  0.878906  3.320312 -1.464844 -1.992187  0.878906
> [6,]  0.000000  0.000000  0.000000  0.000000 -1.464844  2.441406  0.878906 -1.464844
> [7,]  0.000000  0.878906  0.000000 -1.464844 -1.992188  0.878906  3.320313 -1.464844
> [8,]  0.000000  0.000000  0.000000  0.000000  0.878906 -1.464844 -1.464844  2.441406
>
>
> However, this is what I get after switching the kronecker(A,B) to (B,A)
>
> r2          [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]      [,8]
> [1,]  2.441406 -1.464844 -1.464844  0.878906  0.000000  0.000000  0.000000  0.000000
> [2,] -1.464844  2.441406  0.878906 -1.464844  0.000000  0.000000  0.000000  0.000000
> [3,] -1.464844  0.878906  3.320312 -1.992187 -1.464844  0.878906  0.000000  0.000000
> [4,]  0.878906 -1.464844 -1.992187  3.320312  0.878906 -1.464844  0.000000  0.000000
> [5,]  0.000000  0.000000 -1.464844  0.878906  3.320312 -1.992187 -1.464844  0.878906
> [6,]  0.000000  0.000000  0.878906 -1.464844 -1.992188  3.320313  0.878906 -1.464844
> [7,]  0.000000  0.000000  0.000000  0.000000 -1.464844  0.878906  2.441406 -1.464844
> [8,]  0.000000  0.000000  0.000000  0.000000  0.878906 -1.464844 -1.464844  2.441406
>
>
> It is very close. they differ in r2[2,2], r2[3,3] , [7,7] , [6,6] because
> they have been exchanged
> How can we fix this . We are almost there !
> Thanks.
>
> Laz
>
>
>
>
>
> On 6/13/2014 12:24 AM, Kevin Wright wrote:
>
>  I don't have R with me, but what happens if you switch kronecker(A,B)
> with kronecker(B,A)  ?
>
>  kw
>
>
>
> On Thu, Jun 12, 2014 at 10:22 PM, Laz <lmramba@ufl.edu> wrote:
>
>>  Dear Kevin,
>>
>> I have tried your function below and I get slightly different results
>> from my old function. I don't know which one is correct or if I am making a
>> mistake somewhere.
>>
>> Consider a field experiment with a total of 4 genotypes per block. There
>> are 2 blocks each with 2 rows and 2 columns. The total rows=2 and total
>> columns =4. Thus, N=total number of observations=2*4=8. Assume a spatial
>> correlation of 0.3 on both the rows and columns and residual error=1.
>>
>>  b=2;g=4;rb=2;cb=2;r=2;c=4
>>
>>    x y block genotypes
>> 1 1 1     1         1
>> 2 2 1     1         2
>> 3 1 2     1         4
>> 4 2 2     1         3
>> 5 3 1     2         3
>> 6 4 1     2         1
>> 7 3 2     2         2
>> 8 4 2     2         4
>>
>>
>>
>> ### your function
>> #################
>> rspat2<-function(rhox, rhoy, r, c,s2e = 1){
>> require(mvtnorm)
>> sigAR <- diag(r)
>> sigAR <- rhox^abs(row(sigAR) - col(sigAR))
>> sigAC <- diag(c)
>> sigAC <- rhoy^abs(row(sigAC) - col(sigAC))
>> sig <- s2e * kronecker(sigAR, sigAC)
>> round(solve(sig),5)
>> }
>>
>> rspat2(rhox=0.3, rhoy=0.3, r=2,c=4,s2e = 1)
>>
>>
>>  sig       [,1]  [,2]  [,3]   [,4]   [,5]  [,6]  [,7]   [,8]
>> [1,] 1.0000 0.300 0.090 0.0270 0.3000 0.090 0.027 0.0081
>> [2,] 0.3000 1.000 0.300 0.0900 0.0900 0.300 0.090 0.0270
>> [3,] 0.0900 0.300 1.000 0.3000 0.0270 0.090 0.300 0.0900
>> [4,] 0.0270 0.090 0.300 1.0000 0.0081 0.027 0.090 0.3000
>> [5,] 0.3000 0.090 0.027 0.0081 1.0000 0.300 0.090 0.0270
>> [6,] 0.0900 0.300 0.090 0.0270 0.3000 1.000 0.300 0.0900
>> [7,] 0.0270 0.090 0.300 0.0900 0.0900 0.300 1.000 0.3000
>> [8,] 0.0081 0.027 0.090 0.3000 0.0270 0.090 0.300 1.0000
>>
>>  solve(sig)
>>
>>          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
>> [1,]  1.20758 -0.36228  0.00000  0.00000 -0.36228  0.10868  0.00000  0.00000
>> [2,] -0.36228  1.31627 -0.36228  0.00000  0.10868 -0.39488  0.10868  0.00000
>> [3,]  0.00000 -0.36228  1.31627 -0.36228  0.00000  0.10868 -0.39488  0.10868
>> [4,]  0.00000  0.00000 -0.36228  1.20758  0.00000  0.00000  0.10868 -0.36228
>> [5,] -0.36228  0.10868  0.00000  0.00000  1.20758 -0.36228  0.00000  0.00000
>> [6,]  0.10868 -0.39488  0.10868  0.00000 -0.36228  1.31627 -0.36228  0.00000
>> [7,]  0.00000  0.10868 -0.39488  0.10868  0.00000 -0.36228  1.31627 -0.36228
>> [8,]  0.00000  0.00000  0.10868 -0.36228  0.00000  0.00000 -0.36228  1.20758
>>
>>
>>     >
>>
>> However, using my code I get slightly different results shown below:
>> rspat0<-function(rhox,rhoy,s2e=1)
>> {
>>   R<-s2e*eye(N)
>>    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<-solve(R)
>>   IR
>> }
>>
>>    R       [,1]  [,2]   [,3]  [,4]  [,5]   [,6]  [,7]   [,8]
>> [1,] 1.0000 0.300 0.3000 0.090 0.090 0.0270 0.027 0.0081
>> [2,] 0.3000 1.000 0.0900 0.300 0.300 0.0900 0.090 0.0270
>> [3,] 0.3000 0.090 1.0000 0.300 0.027 0.0081 0.090 0.0270
>> [4,] 0.0900 0.300 0.3000 1.000 0.090 0.0270 0.300 0.0900
>> [5,] 0.0900 0.300 0.0270 0.090 1.000 0.3000 0.300 0.0900
>> [6,] 0.0270 0.090 0.0081 0.027 0.300 1.0000 0.090 0.3000
>> [7,] 0.0270 0.090 0.0900 0.300 0.300 0.0900 1.000 0.3000
>> [8,] 0.0081 0.027 0.0270 0.090 0.090 0.3000 0.300 1.0000
>>
>>
>>     >
>>
>>
>>    round(IR,5)         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
>> [1,]  1.20758 -0.36228 -0.36228  0.10868  0.00000  0.00000  0.00000  0.00000
>> [2,] -0.36228  1.31627  0.10868 -0.39488 -0.36228  0.00000  0.10868  0.00000
>> [3,] -0.36228  0.10868  1.20758 -0.36228  0.00000  0.00000  0.00000  0.00000
>> [4,]  0.10868 -0.39488 -0.36228  1.31627  0.10868  0.00000 -0.36228  0.00000
>> [5,]  0.00000 -0.36228  0.00000  0.10868  1.31627 -0.36228 -0.39488  0.10868
>> [6,]  0.00000  0.00000  0.00000  0.00000 -0.36228  1.20758  0.10868 -0.36228
>> [7,]  0.00000  0.10868  0.00000 -0.36228 -0.39488  0.10868  1.31627 -0.36228
>> [8,]  0.00000  0.00000  0.00000  0.00000  0.10868 -0.36228 -0.36228  1.20758
>>
>>
>>     >
>>
>> Sorry for disturbing you again but I thought you can help me in
>> identifying the expected correct Residual structure as your code is more
>> efficient than mine.
>>
>> Thanks very much.
>>
>> Laz
>>
>>
>> On 6/12/2014 7:45 PM, Kevin Wright wrote:
>>
>> There's a much easier way.  See below.
>>
>> Kevin Wright
>>
>>
>> require(mvtnorm)
>> nr <- 4; nc <- 2
>> rho.r <- .1; rho.c <- 0.1
>> sigAR <- diag(nr)
>> sigAR <- rho.r^ abs(row(sigAR) - col(sigAR))
>> sigAC <- diag(nc)
>> sigAC <- rho.c ^ abs(row(sigAC) - col(sigAC))
>> sig <- 1 * kronecker(sigAR, sigAC)
>> round(solve(sig),5)
>>
>> ##          [,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.02030  0.01020 -0.10203  0.00000  0.00000  0.00000
>> 0.00000
>> ## [3,] -0.10203  0.01020  1.03051 -0.10305 -0.10203  0.01020  0.00000
>> 0.00000
>> ## [4,]  0.01020 -0.10203 -0.10305  1.03051  0.01020 -0.10203  0.00000
>> 0.00000
>> ## [5,]  0.00000  0.00000 -0.10203  0.01020  1.03051 -0.10305 -0.10203
>> 0.01020
>> ## [6,]  0.00000  0.00000  0.01020 -0.10203 -0.10305  1.03051  0.01020
>> -0.10203
>> ## [7,]  0.00000  0.00000  0.00000  0.00000 -0.10203  0.01020  1.02030
>> -0.10203
>> ## [8,]  0.00000  0.00000  0.00000  0.00000  0.01020 -0.10203 -0.10203
>> 1.02030
>>
>>
>>
>>
>>  On Thu, Jun 12, 2014 at 5:35 AM, Laz <lmramba@ufl.edu> wrote:
>>
>>> 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.
>>>
>>
>>
>>
>> --
>> Kevin Wright
>>
>>
>>
>
>
> --
> Kevin Wright
>
>
>


-- 
Kevin Wright

	[[alternative HTML version deleted]]

