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