Laz
2014-Jun-12 10:35 UTC
[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]]
ONKELINX, Thierry
2014-Jun-12 11:53 UTC
[R] help in writing an R-function for Residual correlated structures
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 at 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 at r-project.org [mailto:r-help-bounces at r-project.org]
Namens Laz
Verzonden: donderdag 12 juni 2014 12:35
Aan: r-help at 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 at 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.