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


	[[alternative HTML version deleted]]

