x/y/row/col are just arbitrary labels. Change nr to nx. Call them whatever you want. Kevin On Fri, Jun 13, 2014 at 10:16 AM, laz wrote: > The results are similar now but note that in this design, there are > total of 4 columns because x ranges from 1:4 and a total of 2 rows since y > ranges from 1 to 2. > > des1 x y block gen > 1 1 1 1 4 > 3 1 2 1 1 > 2 2 1 1 3 > 4 2 2 1 2 > 5 3 1 2 4 > 7 3 2 2 3 > 6 4 1 2 2 > 8 4 2 2 1 > > > Thus, if we specify correct that nr <- 2; nc <- 4 we still have some > values interchanged whereas if we force nr <- 4; nc <- 2 then the results > are all similar. > How should we fix this one? > > Thanks again Kevin. > > Regards, > Laz > > > > On 6/13/2014 10:52 AM, Kevin Wright wrote: > > 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 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 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 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 > > > -- Kevin Wright [[alternative HTML version deleted]]