Dear R users,
Dear R users,
(I had not included two more functions in the previous mail. This
version is complete)
There is a small problem which I don't know how to sort it out, based on
the former example I had explained earlier own.
I am calling my own functions which are based on simulations as below:
library(gmp)
library(knitr) # load this packages for publishing results
library(matlab)
library(Matrix)
library(psych)
library(foreach)
library(epicalc)
library(ggplot2)
library(xtable)
library(gdata)
library(gplots)
####################################
# function to calculate heritability
herit<-function(varG,varR=1)
{
h<-4*varG/(varG+varR)
h
}
h<-herit(0.081,1);h
###################################
# function to calculate random error
varR<-function(varG,h2)
{
varR<- varG*(4-h2)/h2
varR
}
#system.time(h<-varR(0.081,0.3));h
##########################################
# function to calculate treatment variance
varG<-function(varR=1,h2)
{
varG<-varR*h2/(4-h2)
varG
}
# system.time(h<-varG(1,0.3));h
###############################
# calculating R inverse from spatial data
rspat<-function(rhox=0.6,rhoy=0.6)
{
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
}
### a function to generate A sparse matrix from a pedigree
ZGped<-function(ped)
{
ped2<-data.frame(ped)
lenp2<-length(unique(ped2$V1));lenp2 # how many Genotypes in total in
the pedigree =40
ln2<-length(g);ln2#ln2=nrow(matdf)=30
# calculate the new Z
Zped<-model.matrix(~ matdf$genotypes -1)# has order N*t = 180 by 30
dif<-(lenp2-ln2);dif # 40-30=10
#print(c(lenp2,ln2,dif))
zeromatrix<-zeros(nrow(matdf),dif);zeromatrix # 180 by 10
Z<-cbind(zeromatrix,Zped) # Design Matrix for random effect
(Genotypes): 180 by 40
# calculate the new G
M<-matrix(0,lenp2,lenp2) # 40 by 40
for (i in 1:nrow(ped2)) { M[ped2[i, 1], ped2[i, 2]] <- ped2[i, 3] }
G<-s2g*M # Genetic Variance covariance matrix for pedigree 2: 40 by 40
IG<-solve(G)
results<-c(IG, Z)
results
}
#### Three main functions here #####
### function 1: generate a design (dataframe)
setup<-function(b,g,rb,cb,r,c,h2,rhox=0.6,rhoy=0.6,ped="F")
{
# where
# b = number of blocks
# t = number of treatments per block
# rb = number of rows per block
# cb = number of columns per block
# s2g = variance within genotypes
# h2 = heritability
# r = total number of rows for the layout
# c = total number of columns for the layout
### Check points
if(b==" ")
stop(paste(sQuote("block")," cannot be missing"))
if(!is.vector(g) | length(g)<3)
stop(paste(sQuote("treatments")," should be a vector and more than
2"))
if(!is.numeric(b))
stop(paste(sQuote("block"),"is not of class", sQuote("numeric")))
if(length(b)>1)
stop(paste(sQuote("block"),"has to be only 1 numeric value"))
if(!is.whole(b))
stop(paste(sQuote("block"),"has to be an", sQuote("integer")))
## Compatibility checks
if(rb*cb !=length(g))
stop(paste(sQuote("rb x cb")," should be equal to number of
treatment", sQuote("g")))
if(length(g) != rb*cb)
stop(paste(sQuote("the number of treatments"), "is not equal to",
sQuote("rb*cb")))
## Generate the design
g<<-g
genotypes<-times(b) %do% sample(g,length(g))
#genotypes<-rep(g,b)
block<-rep(1:b,each=length(g))
genotypes<-factor(genotypes)
block<-factor(block)
### generate the base design
k<-c/cb # number of blocks on the x-axis
x<<-rep(rep(1:r,each=cb),k) # X-coordinate
#w<-rb
l<-cb
p<-r/rb
m<-l+1
d<-l*b/p
y<<-c(rep(1:l,r),rep(m:d,r)) # Y-coordinate
## compact
matdf<<-data.frame(x,y,block,genotypes)
N<<-nrow(matdf)
mm<-summ(matdf)
ss<-des(matdf)
## Identity matrices
X<<-model.matrix(~block-1)
h2<<-h2;rhox<<-rhox;rhoy<<-rhoy
s2g<<-varG(varR=1,h2)
## calculate G and Z
ifelse(ped == "F",
c(IG<<-(1/s2g)*eye(length(g)),Z<<-model.matrix(~matdf$genotypes-1)),
c(IG<<- ZGped(ped)[[1]],Z<<-ZGped(ped)[[2]]))
## calculate R and IR
s2e<-1
ifelse(rhox==0 | rhoy==0, IR<<-(1/s2e)*eye(N),
IR<<-rspat(rhox=rhox,rhoy=rhoy))
C11<-t(X)%*%IR%*%X
C11inv<-solve(C11)
K<<-IR%*%X%*%C11inv%*%t(X)%*%IR
return(list( matdf= matdf,summary=mm,description=ss))
}
matrix0<-setup(b=4,g=seq(1,4,1),rb=2,cb=2,r=4,c=4,h2=0.3,rhox=0.6,rhoy=0.6,ped="F")[1]$matdf;
matrix0
x y block genotypes
1 1 1 1 1
2 1 2 1 3
3 2 1 1 2
4 2 2 1 4
5 3 1 2 1
6 3 2 2 3
7 4 1 2 4
8 4 2 2 2
9 1 3 3 1
10 1 4 3 2
11 2 3 3 4
12 2 4 3 3
13 3 3 4 1
14 3 4 4 2
15 4 3 4 3
16 4 4 4 4
### function 2
mainF<-function(criteria=c("A","D"))
{
### Variance covariance matrices
temp<-t(Z)%*%IR%*%Z+IG - t(Z)%*%K%*%Z
C22<-solve(temp)
## calculate trace or determinant
traceI<<-sum(diag(C22)) ## A-Optimality
doptimI<<-log(det(C22)) # D-Optimality
if(criteria=="A") return(traceI)
if(criteria=="D") return(doptimI)
else{return(c(traceI,doptimI))}
}
start0<-mainF(criteria="A");start0
[1] 0.1863854
### function 3 : A function that swaps pairs of treatments randomly
swapsimple<-function(matdf,ped="F")
{
matdf<-as.data.frame(matdf)
attach(matdf,warn.conflict=FALSE)
b1<-sample(matdf$block,1,replace=TRUE);b1
gg1<-matdf$genotypes[block==b1];gg1
g1<-sample(gg1,2);g1
samp<-Matrix(c(g1=g1,block=b1),nrow=1,ncol=3,
dimnames=list(NULL,c("gen1","gen2","block")));samp
newGen<-matdf$genotypes
newG<-ifelse(matdf$genotypes==samp[,1] &
block==samp[,3],samp[,2],matdf$genotypes)
NewG<-ifelse(matdf$genotypes==samp[,2] & block==samp[,3],samp[,1],newG)
NewG<-factor(NewG)
## now, new design after swapping is
newmatdf<-cbind(matdf,NewG)
newmatdf<-as.data.frame(newmatdf)
mm<-summ(newmatdf)
ss<-des(newmatdf)
## Identity matrices
#X<<-model.matrix(~block-1)
#s2g<<-varG(varR=1,h2)
## calculate G and Z
ifelse(ped == "F",
c(IG<<-(1/s2g)*eye(length(g)),Z<<-model.matrix(~newmatdf$NewG-1)),
c(IG<<- ZGped(ped)[[1]],Z<<-ZGped(ped)[[2]]))
## calculate R and IR
C11<-t(X)%*%IR%*%X
C11inv<-solve(C11)
K<<-IR%*%X%*%C11inv%*%t(X)%*%IR
#Nmatdf<-newmatdf[,c(1,2,3,5)]
names(newmatdf)[names(newmatdf)=="genotypes"] <- "old_G"
names(newmatdf)[names(newmatdf)=="NewG"] <- "genotypes"
#newmatdf <- remove.vars(newmatdf, "old_G")
newmatdf$old_G <- newmatdf$old_G <- NULL
#matdf<-newmatdf
newmatdf
}
matdf<-swapsimple(matdf,ped="F")
matdf
x y block genotypes
1 1 1 1 1
2 1 2 1 3
3 2 1 1 2
4 2 2 1 4
5 3 1 2 4
6 3 2 2 3
7 4 1 2 1
8 4 2 2 2
9 1 3 3 1
10 1 4 3 2
11 2 3 3 4
12 2 4 3 3
13 3 3 4 1
14 3 4 4 2
15 4 3 4 3
16 4 4 4 4
which(matrix0$genotypes != matdf$genotypes)
[1] 5 7
# This is fine because I expected a maximum of 1 pair to change, so I
have a maximum of 2 positions swapped on the first iteration.
# If I swap 10 times (iterations=10), I expect a maximum of 20
positions to change
### The final function (where I need your help more )
fun <- function(n = 10){
matrix0<-setup(b=4,g=seq(1,4,1),rb=2,cb=2,r=4,c=4,h2=0.3,rhox=0.6,rhoy=0.6,ped="F")[1]$matdf
# matrix0 is the original design before swapping any pairs of treatments
res <- list(mat = NULL, Design_best = matrix0, Original_design =
matrix0)
start0<-mainF(criteria="A")
# start0 is the original trace
res$mat <- rbind(res$mat, c(trace = start0, iterations = 0))
for(i in seq_len(n)){
# now swap the pairs of treatments from the original design, n times
matdf<-swapsimple(matdf,ped="F")
if(mainF(criteria="A") < start0){
start0<- mainF(criteria="A")
res$mat <- rbind(res$mat, c(trace = start0, iterations = i))
res$Design_best <- matdf
}
}
res
}
res<-fun(50)
res
$mat
trace iterations
[1,] 0.1938285 0
[2,] 0.1881868 1
[3,] 0.1871099 17
[4,] 0.1837258 18
[5,] 0.1812291 19
### here is the problem
which(res$Design_best$genotypes != res$Original_design$genotypes) #
always get a pair of difference
[1] 2 3 4 5 6 7 10 11 13 14 15 16
## I expect a maximum of 8 changes but I get 12 changes which means that
function only dropped the traces when trace_j > trace_i but did not drop
the design !!
How do I fix this ?????
Kind regards,
lazarus
On 10/19/2013 5:03 PM, Rui Barradas wrote:
Hello,
Seems simple.
fun <- function(n = 10){
matd <- matrix(sample(1:30,30, replace=FALSE), ncol=5, nrow=6)
res <- list(mat = NULL, Design_best = matd, Original_design = matd)
trace <- sum(diag(matd))
res$mat <- rbind(res$mat, c(trace = trace, iterations = 0))
for(i in seq_len(n)){
matd <- matrix(sample(1:30,30, replace=FALSE), ncol=5, nrow=6)
if(sum(diag(matd)) < trace){
trace <- sum(diag(matd))
res$mat <- rbind(res$mat, c(trace = trace, iterations = i))
res$Design_best <- matd
}
}
res
}
fun()
fun(20)
Hope this helps,
Rui Barradas
Em 19-10-2013 18:41, laz escreveu:
Dear R users,
Suppose I want to randomly generate some data, in matrix form, randomly
swap some of the elements and calculate trace of the matrix for each of
these stages. If the value of trace obtained in the later is bigger than
the former, drop the latter matrix and go back to the former matrix,
swap some elements of the matrix again and calculate the trace. If the
recent trace is smaller than the previous one, accept the matrix as the
current . Use the current matrix and swap elements again. repeat the
whole process for a number of times, say, 10. The output from the
function should display only the original matrix and its value of trace,
trace values of successful swaps and their iteration counts and the
final best matrix that had the smallest value of trace, together with
its trace value.
For example
## original
> matd<-matrix(sample(1:30,30,replace=FALSE),ncol=5,nrow=6,byrow=FALSE)
> matd
[,1] [,2] [,3] [,4] [,5]
[1,] 12 27 29 16 19
[2,] 25 10 7 22 13
[3,] 14 23 3 11 21
[4,] 28 6 5 2 18
[5,] 24 20 1 17 26
[6,] 9 4 30 8 15
> trace<-sum(diag(matd))
> trace
[1] 53
# 1st iteration
[,1] [,2] [,3] [,4] [,5]
[1,] 24 29 20 25 17
[2,] 16 1 30 9 5
[3,] 18 22 2 10 26
[4,] 23 27 19 21 28
[5,] 15 6 8 3 13
[6,] 12 14 7 11 4
> trace<-sum(diag(matd))
> trace
[1] 61
## drop this matrix because 61 > 53
# 2nd iteration
> matd<-matrix(sample(1:30,30,replace=FALSE),ncol=5,nrow=6,byrow=FALSE)
> matd
[,1] [,2] [,3] [,4] [,5]
[1,] 2 28 23 15 14
[2,] 27 9 10 29 7
[3,] 5 18 12 1 11
[4,] 8 4 30 16 24
[5,] 25 19 26 6 13
[6,] 17 22 3 20 21
> trace<-sum(diag(matd))
> trace
[1] 52
## accept this matrix because 52 < 53
### 3rd iteration
> matd<-matrix(sample(1:30,30,replace=FALSE),ncol=5,nrow=6,byrow=FALSE)
> matd
[,1] [,2] [,3] [,4] [,5]
[1,] 1 29 17 8 6
[2,] 21 23 10 7 14
[3,] 22 4 12 26 9
[4,] 3 13 11 30 15
[5,] 5 24 18 16 2
[6,] 20 25 19 27 28
> trace<-sum(diag(matd))
> trace
[1] 68
## drop this matrix because 68 > 52
## 4th iteration
> matd<-matrix(sample(1:30,30,replace=FALSE),ncol=5,nrow=6,byrow=FALSE)
> matd
[,1] [,2] [,3] [,4] [,5]
[1,] 2 6 5 28 15
[2,] 9 12 13 19 24
[3,] 3 22 14 11 29
[4,] 30 20 17 7 23
[5,] 18 27 21 1 10
[6,] 25 16 4 8 26
> trace<-sum(diag(matd))
> trace
[1] 45
## accept this matrix because 45 < 52
The final results will be:
$mat
trace iterations
[1,] 53 0
[2,] 52 2
[3,] 45 4
$ Design_best
[,1] [,2] [,3] [,4] [,5]
[1,] 2 6 5 28 15
[2,] 9 12 13 19 24
[3,] 3 22 14 11 29
[4,] 30 20 17 7 23
[5,] 18 27 21 1 10
[6,] 25 16 4 8 26
$ Original_design
[,1] [,2] [,3] [,4] [,5]
[1,] 12 27 29 16 19
[2,] 25 10 7 22 13
[3,] 14 23 3 11 21
[4,] 28 6 5 2 18
[5,] 24 20 1 17 26
[6,] 9 4 30 8 15
Regards,
Laz
______________________________________________
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.
