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
[email protected]
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: [email protected] [mailto:[email protected]] Namens
Laz
Verzonden: donderdag 12 juni 2014 12:35
Aan: [email protected]
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]]
______________________________________________
[email protected] 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.
______________________________________________
[email protected] 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.
