[R] loop not working my way
Rui Barradas
ruipbarradas at sapo.pt
Sun Oct 20 18:54:01 CEST 2013
Hello,
Make a _simple_ example, I don't see what packages like knitr or ggplot2
have anything to do with your problem.
Like this is, I think you're asking too much from r-help.
Rui Barradas
Em 19-10-2013 23:38, Laz escreveu:
> 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 at 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.
>>
>
>
More information about the R-help
mailing list