Dear R users,
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 three 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 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:09 PM, Laz wrote:
> Thank you so very much!
> It works like a charm !!!
>
> Regards,
> Laz
>
> On 10/19/2013 5:03 PM, Rui Barradas wrote:
>> 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)
>
[[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.
