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
>>
>> ______________________________________________
>> [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.
>
[[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.
