R Under development (unstable) (2023-11-12 r85514 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # plotMonitor <- function(x,...){
> #   mmin <- min(x)
> #   mmax <- max(x)
> #   plot(x[1,], type="l", ylim=c(mmin,mmax))
> #   if(nrow(x) > 1){
> #     for(i in 2:nrow(x)){
> #       par(new=TRUE)
> #       plot(x[i,], col=i, ylim=c(mmin,mmax), ylab="",xlab="", type="l",...)
> #     }
> #   }
> #   legend("topright",legend = rownames(x), col=1:(nrow(x)), bty="n", lty=1)
> # }
> # 
> # library(sommer)
> # data("DT_cpdata")
> # DT <- DT_cpdata
> # GT <- GT_cpdata
> # MP <- MP_cpdata
> # DT$Yield <- imputev(DT$Yield)
> # #### create the variance-covariance matrix
> # K <- A.mat(GT) # additive relationship matrix
> # K <- K + diag(1e-3,nrow(K))
> # #### look at the data and fit the model
> # head(DT)
> # mix1 <- mmer(Yield~1,
> #              random=~vsr(id,Gu=K)
> #              + Rowf,
> #              rcov=~units,
> #              data=DT)
> # summary(mix1)$varcomp
> # 
> # Ki <- as(solve(K), Class="sparseMatrix")
> # mix2 <- mmec(fixed=Yield~1,
> #              random=~vsc(ids(id),Gu=Ki) + Rowf,
> #              rcov=~units, return.param = F,
> #              data=DT)
> # 
> # mix2$monitor[,15]
> # plotMonitor(mix2$monitor)
> # pp <- predict.mmec(object=mix2,classify = c("id","(Intercept)"))#, ignore="Rowf")
> # head(pp$pvals)
> # 
> # Z <- list(model.matrix(~id-1, data=DT),model.matrix(~Rowf-1,data=DT)
> # )
> # Zind <- 1:2
> # 
> # A <- list(
> #   K, 
> #   diag(13)
> # ) #
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # thetaI <- list(
> #   matrix(1000,1,1),
> #   matrix(1000,1,1),
> #   matrix(1000,1,1)
> # );thetaI
> # 
> # thetaC <- list(
> #   matrix(1,1,1),
> #   matrix(1,1,1),
> #   matrix(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~1, data=DT)
> # 
> # y <- as.matrix(DT$Yield)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # 
> # ## apply the function
> # weightInf=rep(1,40) # weights for the information matrix
> # weightEmInf=c(seq(.9,.1,-.2),rep(0,36));weightEmInf # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # # compare results
> # res3$monitor
> # summary(mix1)$varcomp
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # # CS + DIAGONAL MODEL
> # 
> # library(sommer)
> # data("DT_example")
> # DT <- DT_example
> # K <- A_example
> # #### look at the data and fit the model
> # head(DT)
> # DT$Yield <- scale(DT$Yield)
> # mix1 <- mmer(Yield~Env,
> #              random= ~Name + vsr(dsr(Env),Name),
> #              rcov= ~ vsr(dsr(Env),units),
> #              data=DT)
> # summary(mix1)$varcomp
> # 
> # mix2 <- mmec(Yield~Env,
> #              random= ~ Name + vsc(dsr(Env),ids(Name)),
> #              rcov= ~ vsc(dsr(Env),ids(units)),
> #              return.param = F,
> #              data=DT)
> # plotMonitor(mix2$monitor)
> # summary(mix1)$varcomp
> # mix2$sigma
> # 
> # zz <- with(DT, vsr(dsr(Env),Name))
> # 
> # Z <- c(list(model.matrix(~Name-1, data=DT)),zz$Z)
> # 
> # Zind <- c(1,2,2,2)
> # 
> # A <- list(diag(41), diag(41))#rep(list(diag(41)),4)
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # thetaI <- list(
> #   matrix(10,1,1),
> #   diag(10,3,3),
> #   diag(10,3,3)
> # );thetaI
> # 
> # thetaC <- list(
> #   matrix(1,1,1),
> #   diag(1,3,3),
> #   diag(1,3,3)
> # );thetaC
> # 
> # X <- model.matrix(~Env, data=DT)
> # 
> # y <- as.matrix(DT$Yield)
> # 
> # DTx <- DT; DTx$units <- as.factor(1:nrow(DTx))
> # ss <- with(DTx, vsr(dsr(Env),units) )
> # 
> # S <- ss$Z #list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # 
> # ## apply the function
> # weightInf=rep(1,40); # weights for the information matrix
> # weightEmInf=c(seq(.9,.1,-.1),rep(0,36)); # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # # compare results
> # res3$monitor
> # summary(mix1)$varcomp
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # # UNSTRUCTURED MODEL
> # library(sommer)
> # data("DT_example")
> # DT <- DT_example
> # K <- A_example
> # #### look at the data and fit the model
> # head(DT)
> # DT$Yield <- scale(DT$Yield)
> # mix1 <- mmer(Yield~Env,
> #              random= ~ vsr(usr(Env),Name),
> #              rcov= ~ vsr(dsr(Env),units),
> #              data=DT)
> # summary(mix1)$varcomp[,1]
> # mix1$Beta
> # 
> # # > summary(mix1)$varcomp[,1]
> # # [1] 0.76660279 0.29901471 0.22166682 0.31243923 0.01923234 0.42072739 0.24321083 0.27761802
> # # [9] 0.12513727
> # 
> # mix2 <- mmec(fixed=Yield~Env-1,
> #               random= ~ vsc(usr(Env),ids(Name)),
> #               rcov= ~ vsc(dsr(Env),ids(units)),
> #               return.param = F,
> #              # emweight = rep(1,30),
> #              nIters = 30,
> #               data=DT)
> # mix2$theta
> # plotMonitor(mix2$monitor, cex=1)
> # mix2$sigma
> # mix2$monitor[,12]
> # plot(summary(mix1)$varcomp[,1],mix2$sigma)
> # 
> # zz <- with(DT, vsr(dsr(Env),Name))
> # 
> # Z <- zz$Z
> # 
> # Zind <- rep(1,length(Z))
> # 
> # A <- list( diag(41) )#rep(list(diag(41)),4)
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # tt = ((unsm(3)/5) + diag(.8,3,3)) * 10 ;# tt[lower.tri(tt)]=0;tt
> # thetaI <- list(
> #   tt,
> #   diag(10,3,3)
> # );thetaI
> # 
> # ttc= unsm(3);ttc[lower.tri(ttc)]=0;ttc
> # thetaC <- list(
> #   ttc,
> #   diag(1,3,3)
> # );thetaC
> # 
> # X <- model.matrix(~Env, data=DT)
> # 
> # y <- as.matrix(DT$Yield)
> # 
> # DTx <- DT; DTx$units <- as.factor(1:nrow(DTx))
> # ss <- with(DTx, vsr(dsr(Env),units) )
> # 
> # S <- ss$Z #list(diag(length(y)))
> # 
> # Sind <- rep(1,3)
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # 
> # ## apply the function
> # weightInf=rep(1,50);weightInf # weights for the information matrix
> # weightEmInf=rep(1,50)
> # weightEmInf=c(seq(1,.1,-.1),rep(0,50));weightEmInf # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # summary(mix1)$varcomp
> # res3$monitor
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # 
> # 
> # library(sommer)
> # data(DT_technow)
> # DT <- DT_technow
> # Md <- Md_technow
> # Mf <- Mf_technow
> # Md <- (Md*2) - 1
> # Mf <- (Mf*2) - 1
> # Ad <- A.mat(Md)
> # Af <- A.mat(Mf)
> # DT$GY <- scale(DT$GY)
> # ####=========================================####
> # ####=========================================####
> # mix1 <- mmer(GY~1,
> #              random=~vsr(dent,Gu=Ad) + vsr(flint,Gu=Af),
> #              rcov=~units,
> #              data=DT)
> # summary(mix1)$varcomp[,1]
> # 
> # Adi <- as(solve(Ad + diag(ncol(Ad))*1e-6 ),Class="sparseMatrix")
> # Afi <- as(solve(Af + diag(ncol(Af))*1e-6 ),Class="sparseMatrix")
> # 
> # mix2 <- mmec(GY~1,
> #              random=~vsc(ids(dent),Gu=Adi) + vsc(ids(flint),Gu=Afi),
> #              rcov=~units, return.param = F,
> #              data=DT)
> # mix2$monitor
> # mix2$sigma
> # 
> # z1=model.matrix(~dent-1, data=DT); colnames(z1) <- gsub("dent","",colnames(z1))
> # z2=model.matrix(~flint-1, data=DT); colnames(z2) <- gsub("flint","",colnames(z2))
> # 
> # Z <- list(z1,z2)
> # 
> # Zind <- 1:2
> # 
> # A <- list(
> #   Ad[colnames(z1),colnames(z1)], 
> #   Af[colnames(z2),colnames(z2)]
> # ) 
> # A <- lapply(A, function(x){x + diag(1e-3,ncol(x),ncol(x))})
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # thetaI <- list(
> #   matrix(10,1,1),
> #   matrix(10,1,1),
> #   matrix(10,1,1)
> # );thetaI
> # 
> # thetaC <- list(
> #   matrix(1,1,1),
> #   matrix(1,1,1),
> #   matrix(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~1, data=DT)
> # 
> # y <- as.matrix(DT$GY)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # ## apply the function
> # weightInf=rep(1,40);weightInf # weights for the information matrix
> # weightEmInf=c(seq(.9,.1,-.2),rep(0,36));weightEmInf # weights for the EM information matrix
> # # v2=rep(1,40)
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=FALSE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # summary(mix1)$varcomp
> # res3$monitor
> # 
> # dim(res3$Ci)
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # 
> # library(sommer)
> # data(DT_ige)
> # DT <- DT_ige
> # Af <- A_ige
> # An <- A_ige
> # DT$trait <- scale(imputev(DT$trait))
> # ### Direct genetic effects model
> # mix1 <- mmer(trait ~ block,
> #              random = ~vsr(focal,Gu=Af) + vsr(neighbour,Gu=An),
> #              rcov = ~ units,
> #              data = DT)
> # summary(mix1)$varcomp
> # 
> # Afi <- as(solve(Af + diag(ncol(Af))*1e-6 ),Class="sparseMatrix")
> # Ani <- as(solve(An + diag(ncol(An))*1e-6 ),Class="sparseMatrix")
> # 
> # mix2 <- mmec(trait ~ block,
> #              random = ~vsc(ids(focal),Gu=Afi) + vsc(ids(neighbour),Gu=Ani),
> #              rcov = ~ units, return.param = F,
> #              data = DT)
> # mix2$monitor
> # plotMonitor(mix2$monitor)
> # 
> # g=(unsm(2)*-.15) + (diag(2)*.4)
> # g
> # mix3 <- mmec(trait ~ block,
> #              random = ~vsc(ids(focal),ids(neighbour),Gu=Afi, meN=2, meThetaC = unsm(2)),# + vsc(ids(neighbour),Gu=Ani),
> #              rcov = ~ units, return.param = F, emweight = rep(1,20),
> #              tolParConv = 1e-6, nIters=40,
> #              data = DT)
> # mix3$monitor
> # mix3$sigma
> # mix3$rTermsNames
> # plotMonitor(mix3$monitor)
> # 
> # 
> # z1=model.matrix(~focal-1,data=DT); colnames(z1) <- gsub("focal","",colnames(z1))
> # z2=model.matrix(~neighbour-1,data=DT);  colnames(z2) <- gsub("neighbour","",colnames(z2))
> # 
> # Z <- list(z1,z2)
> # 
> # Zind <- 1:2
> # 
> # A <- list( 
> #   Af[colnames(z1),colnames(z1)],
> #   An[colnames(z2),colnames(z2)]
> # )
> # 
> # A <- lapply(A, function(x){x + diag(1e-3,ncol(x),ncol(x))})
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # tt = diag(1)*10000
> # thetaI <- list(
> #   tt,
> #   tt,
> #   diag(10000,1,1)
> # );thetaI
> # 
> # ttc=diag(1)
> # thetaC <- list(
> #   ttc,
> #   ttc,
> #   diag(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~block, data=DT)
> # 
> # y <- as.matrix(DT$trait)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # 
> # ## apply the function
> # weightInf=rep(1,40);weightInf # weights for the information matrix
> # weightEmInf=c(seq(.9,.1,-.1),rep(0,36));weightEmInf # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # # compare results
> # summary(mix1)$varcomp
> # res3$monitor
> # 
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # # UNSTRUCTURED MODEL
> # 
> # library(sommer)
> # data(DT_ige)
> # DT <- DT_ige
> # Af <- A_ige
> # An <- A_ige
> # DT$trait <- imputev(DT$trait)
> # ### Direct genetic effects model
> # xx <- c(rep(1,4),rep(0,40))
> # mix1 <- mmer(trait ~ block,
> #              random = ~ gvs(focal, neighbour, Gu=list(Af,Af)), 
> #              rcov = ~ units, iters=10,tolpar = 1e-4,
> #              data = DT, emupdate = xx)
> # summary(mix1)$varcomp
> # 
> # z1=model.matrix(~focal-1,data=DT); colnames(z1) <- gsub("focal","",colnames(z1))
> # z2=model.matrix(~neighbour-1,data=DT);  colnames(z2) <- gsub("neighbour","",colnames(z2))
> # 
> # Z <- list(z1,z2)
> # 
> # Zind <- rep(1,2)
> # A <- list( 
> #   Af
> # )
> # 
> # A <- lapply(A, function(x){x + diag(1e-3,ncol(x),ncol(x))})
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # tt = ((-unsm(2)/5) + diag(.8,2,2)) * 1000 ; #tt[lower.tri(tt)]=0;tt
> # # tt = diag(2)*10000
> # thetaI <- list(
> #   tt,
> #   diag(10000,1,1)
> # );thetaI
> # 
> # ttc= unsm(2);ttc[lower.tri(ttc)]=0;ttc
> # # ttc=diag(2)
> # thetaC <- list(
> #   ttc,
> #   diag(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~block, data=DT)
> # 
> # y <- as.matrix(DT$trait)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # 
> # ## apply the function
> # 
> # weightInf=c(rep(1,100));weightInf
> # weightEmInf=c(rep(1,4),seq(1,.1,-.2),rep(0.05,3),rep(0.05,30));weightEmInf # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # # compare results
> # summary(mix1)$varcomp
> # res3$monitor
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # # UNSTRUCTURED MODEL
> # 
> # library(orthopolynom)
> # data(DT_legendre)
> # DT <- DT_legendre
> # DT$Y <- scale(imputev(DT$Y))
> # mix1<-mmer(Y~ 1 + Xf
> #            , random=~ vsr(usr(leg(X,1)),SUBJECT)
> #            , rcov=~vsr(units), tolpar = 1e-6
> #            , data=DT)
> # summary(mix1)$varcomp
> # 
> # mix2<-mmec(Y~ 1 + Xf
> #            , random=~ vsc(usr(leg(X,1)),ids(SUBJECT))
> #            , rcov=~units,return.param = F
> #            , data=DT)
> # mix2$monitor
> # mix2$sigma
> # 
> # xx=with(DT,vsr(usr(leg(X,1)),SUBJECT))
> # 
> # Z <- xx$Z[c(1,3)]
> # 
> # Zind <- rep(1,2)
> # 
> # A <- list(
> #   diag(100)
> # ) #
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # tt = ((unsm(2)/5) + diag(.8,2,2))  ; #tt[lower.tri(tt)]=0;tt
> # thetaI <- list(
> #   tt,
> #   diag(1,1,1)
> # );thetaI
> # 
> # ttc= unsm(2);ttc[lower.tri(ttc)]=0;ttc
> # # ttc=diag(2)
> # thetaC <- list(
> #   ttc,
> #   diag(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~1 + Xf, data=DT)
> # 
> # y <- as.matrix(DT$Y)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # 
> # ## apply the function
> # weightInf=rep(1,40) # weights for the information matrix
> # weightEmInf=rep(1,40)
> # weightEmInf=c(seq(1,.1,-.2),rep(0,36));weightEmInf # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # # compare results
> # summary(mix1)$varcomp
> # res3$monitor
> # 
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # 
> # data(DT_cornhybrids)
> # DT <- DT_cornhybrids
> # DT<-DT[which(!is.na(DT$Yield)),]
> # nrow(DT)
> # DTi <- DTi_cornhybrids
> # GT <- GT_cornhybrids
> # hybrid2 <- DT # extract cross data
> # A <- GT
> # K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)
> # K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)
> # S <- kronecker(K1, K2) ; dim(S)
> # rownames(S) <- colnames(S) <- levels(hybrid2$SCA)
> # 
> # hybrid2$Yield <- scale(hybrid2$Yield)
> # 
> # mix1 <- mmer(Yield ~ Location,
> #              random = ~ vsr(GCA1,Gu=K1) + vsr(GCA2,Gu=K2),# + vsr(SCA,Gu=S),
> #              rcov=~units,
> #              data=hybrid2)
> # summary(mix1)$varcomp
> # 
> # K1i <- as(solve(K1 + diag(ncol(K1))*1e-6 ),Class="sparseMatrix")
> # K2i <- as(solve(K2 + diag(ncol(K2))*1e-6 ),Class="sparseMatrix")
> # Si <- as(solve(S + diag(ncol(S))*1e-6 ),Class="sparseMatrix")
> # 
> # mix2 <- mmec(Yield ~ Location,
> #              random = ~ vsc(ids(GCA1),Gu=K1i) + vsc(ids(GCA2),Gu=K2i),# + vsX(ids(SCA),Gu=Si),
> #              rcov=~units, return.param = F,
> #              data=hybrid2)
> # mix2$monitor
> # mix2$sigma
> # 
> # z1 <- model.matrix(~GCA1-1,data=DT);colnames(z1) <- gsub("GCA1","",colnames(z1))
> # z2 <- model.matrix(~GCA2-1,data=DT);colnames(z2) <- gsub("GCA2","",colnames(z2))
> # z3 <- model.matrix(~SCA-1,data=DT);colnames(z3) <- gsub("SCA","",colnames(z3))
> # 
> # Z <- list(
> #   z1,z2,z3
> # )
> # 
> # Zind <- 1:3
> # 
> # A <- list(
> #   K1[colnames(z1),colnames(z1)], K2[colnames(z2),colnames(z2)], S[colnames(z3),colnames(z3)]
> # ) #
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # thetaI <- list(
> #   diag(1)*10,
> #   diag(1)*10,
> #   diag(1)*10,
> #   matrix(50,1,1)
> # );thetaI
> # 
> # thetaC <- list(
> #   diag(1),
> #   diag(1),
> #   diag(1),
> #   matrix(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~Location, data=DT)
> # 
> # y <- as.matrix(DT$Yield)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # 
> # ## apply the function
> # weightInf=rep(1,40) # weights for the information matrix
> # weightEmInf=rep(1,40)
> # weightEmInf=c(seq(.9,.5,-.1),rep(0,36));weightEmInf # weights for the EM information matrix
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=TRUE,
> #                       nIters=20, tolParConv=1e-5,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # # compare results
> # summary(mix1)$varcomp
> # res3$monitor
> # 
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # ###############################################################################
> # 
> # library(sommer)
> # 
> # data(DT_h2)
> # DT <- DT_h2
> # 
> # head(DT)
> # length(table(DT$Env))
> # nrow(DT)
> # DT$y <- scale(DT$y)
> # 
> # ans1 <- mmer(y~Env,
> #              random=~vsr(dsr(Env),Name) + vsr(dsr(Env),Block),
> #              rcov=~vsr(dsr(Env),units),
> #              data=DT)
> # summary(ans1)$varcomp
> # 
> # ans2 <- mmec(y~Env,
> #              random=~vsc(dsr(Env),ids(Name)) + vsc(dsr(Env),ids(Block)),
> #              rcov=~vsc(dsr(Env),ids(units)), return.param = F,
> #              # emweight = c(1,rep(0,20)),
> #              # stepweight = c(1,1,1,rep(0,30)),
> #              nIters=30,
> #              data=DT)
> # 
> # plotMonitor(ans2$monitor, cex=.1)
> # ans2$sigma
> # 
> # d=ans2$monitor
> # plotMonitor(d)
> # plot(d[,(ncol(d)-1)], summary(ans1)$varcomp[,1])
> # 
> # DT2 <- DT[with(DT, order(Env)), ]
> # library(asreml)
> # ans3 <- asreml(y~Env,
> #               random=~diag(Env):Name + diag(Env):Block,
> #               residual=~dsum(~units | Env),
> #               
> #               data=DT2)
> # 
> # plot(d[,(ncol(d)-1)], summary(ans3)$varcomp[,1])
> # 
> # 
> # xx=with(DT,vsr(dsr(Env),Name))
> # 
> # Z <- xx$Z
> # 
> # Zind <- rep(1,length(Z))
> # 
> # A <- list(
> #   
> #   diag(41)
> # ) #
> # 
> # Ai <- lapply(A, function(x){solve(x)})
> # 
> # thetaI <- list(
> #   diag(15)*5,#+diag(rnorm(15)),
> #   matrix(10,1,1)
> # );thetaI
> # 
> # thetaC <- list(
> #   diag(15),
> #   matrix(1,1,1)
> # );thetaC
> # 
> # X <- model.matrix(~Env, data=DT)
> # 
> # y <- as.matrix(DT$y)
> # 
> # S <- list(diag(length(y)))
> # 
> # H <- diag(length(y))
> # 
> # addScaleParam <- 0
> # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
> # nn2 <- sum(nn[1:max(Zind)])
> # ff <- diag(nn2)
> # thetaF <- cbind(ff,matrix(0,nn2,1))
> # thetaF
> # 
> # ## apply the function
> # weightInf=rep(1,40) # weights for the information matrix
> # 
> # weightEmInf =c(seq(1,.1,-.1),rep(.01,30))
> # # weightEmInf =rep(.01,30)
> # # weightInfEMv=c(rep(1,8),rep(0,36));weightInfEMv # weights for the EM information matrix
> # 
> # 
> # Z <- lapply(Z, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # S <- lapply(S, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # Ai <- lapply(Ai, function(x){
> #   return(as(x, Class = "sparseMatrix"))
> # })
> # 
> # H <- as(H, Class = "sparseMatrix")
> # 
> # X <- as(X, Class = "sparseMatrix")
> # 
> # y <- as(y, Class = "sparseMatrix")
> # 
> # tt=system.time( 
> #   expr = res3<-ai_mme_sp(X=X,Z=Z, Zind=Zind, 
> #                       Ai=Ai,y=y,
> #                       S=S, 
> #                       H=H, useH=FALSE,
> #                       nIters=15, tolParConv=1e-4,
> #                       tolParInv=1e-6,thetaI=thetaI, 
> #                       thetaC=thetaC,thetaF=thetaF, 
> #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
> #                       weightInf = weightInf,
> #                       verbose=TRUE
> #                       
> #   )
> # )
> # 
> # plot(res3$monitor[,ncol(res3$monitor)], summary(ans1)$varcomp[,1] )
> # 
> # plot(res3$llik[1,])
> # 
> # #######################################################
> # 
> # data(DT_yatesoats)
> # DT <- DT_yatesoats
> # head(DT)
> # DT$Y <- scale(DT$Y)
> # m3 <- mmer(fixed=Y ~ V + N + V:N,
> #            # random = ~ B + B:MP,
> #            rcov=~units,
> #            data = DT)
> # summary(m3)$varcomp
> # 
> # m4 <- mmec(fixed=Y ~ V + N + V:N,
> #            # random = ~ B + B:MP,
> #            rcov=~units, return.param = F,
> #            data = DT)
> # m4$sigma
> # plotMonitor(m4$monitor)
> # str(m4)
> # 
> # 
> # ####################################################
> # 
> # library(sommer)
> # data("DT_cpdata")
> # DT <- DT_cpdata
> # GT <- GT_cpdata
> # MP <- MP_cpdata
> # DT$Yield <- imputev(DT$Yield)
> # #### create the variance-covariance matrix
> # K <- A.mat(GT) # additive relationship matrix
> # K <- K + diag(1e-3,nrow(K))
> # #### look at the data and fit the model
> # head(DT)
> # mix1 <- mmer(Yield~1,
> #              random=~vsr(id,Gu=K)
> #              + Rowf,
> #              rcov=~units,
> #              data=DT)
> # summary(mix1)$varcomp
> # 
> # # Ki <- as(solve(K), Class="sparseMatrix")
> # m <- GT[,1:300]
> # mix2 <- mmec(fixed=Yield~1,
> #               random=~vsc(ids(m)) + Rowf,
> #               rcov=~units, return.param = F,
> #               data=DT)
> # mix2$monitor
> # summary(mix1)$varcomp
> # 
> # 
> # 
> # ?overlay
> # 
> # data("DT_halfdiallel")
> # DT <- DT_halfdiallel
> # head(DT)
> # DT$femalef <- as.factor(DT$female)
> # DT$malef <- as.factor(DT$male)
> # DT$genof <- as.factor(DT$geno)
> # 
> # A <- diag(7); colnames(A) <- rownames(A) <- 1:7;A # if you want to provide a covariance matrix
> # #### model using overlay
> # modh <- mmer(sugar~1, 
> #              random=~vsr(overlay(femalef,malef), Gu=A) 
> #              + genof,
> #              data=DT)
> # 
> # Ai <- as(solve(A), Class="sparseMatrix")
> # modh <- mmec(sugar~1, 
> #              random=~vsc(ids(overlay(femalef,malef)), Gu=Ai) 
> #              + genof, #return.param = T,
> #              data=DT)
> # modh$monitor
> # 
> # #####################################################
> # 
> # ?DT_cpdata
> # 
> # data(DT_cpdata)
> # DT <- DT_cpdata
> # GT <- GT_cpdata
> # MP <- MP_cpdata
> # #### create the variance-covariance matrix
> # A <- A.mat(GT) # additive relationship matrix
> # Ai <- as( solve(A+diag(1e-5,ncol(A),ncol(A))), Class = "sparseMatrix")
> # 
> # head(DT)
> # DT <- DT[,-c(7:8)]
> # DT$color <- as.vector(scale(DT$color))
> # DT$Yield <- as.vector(scale(DT$Yield))
> # head(DT)
> # DT2 <- reshape(DT, idvar = "id", varying = list(5:6),
> #                v.names = "y", direction = "long", timevar = "trait", times =colnames(DT)[5:6] )
> # DT2$trait <- as.factor(DT2$trait)
> # head(DT2)
> # 
> # g=diag(2)*.05 + matrix(.1,2,2);g
> # g2 <- diag(2)*.75;g2
> # mix1 <- mmec(y~trait-1,
> #              random=~vsc(usr(trait, theta = g),ids(id),Gu=Ai),
> #              stepweight = ss,
> #              emweight =rep(1,30),
> #              return.param = F, tolParConv = 1e-6,
> #              rcov=~vsc(dsr(trait),ids(units)),
> #              data=DT2)
> # 
> # mix1$theta
> # mix1$thetaC
> # mix1$sigma
> # mix1$b
> # plot(mix1$u)
> # plot(mix1$llik[1,])
> # 
> # # > summary(mix2)$varcomp
> # # VarComp  VarCompSE    Zratio Constraint
> # # u:id.color-color    0.7557346 0.15068202  5.015426   Positive
> # # u:id.color-Yield    0.0848963 0.07171699  1.183768   Unconstr
> # # u:id.Yield-Yield    0.1394332 0.07012037  1.988484   Positive
> # # u:units.color-color 0.3779553 0.04195655  9.008254   Positive
> # # u:units.Yield-Yield 0.8808832 0.07522796 11.709519   Positive
> # # # 
> # 
> 
> proc.time()
   user  system elapsed 
   0.18    0.04    0.21