Nothing
## ----initialsetup, include=FALSE----------------------------------------------
knitr::opts_chunk$set(cache=FALSE)
## ----eval=FALSE---------------------------------------------------------------
# library(SFSI)
# if(requireNamespace("BGLR")){
# data(wheat, package="BGLR") # Load data from the BGLR package
# }
#
# # Select the environment 1 to work with
# y <- as.vector(scale(wheat.Y[,1]))
#
# # Calculate G matrix
# G <- tcrossprod(scale(wheat.X))/ncol(wheat.X)
#
# # Save data
# save(y, G, file="geno_pheno.RData")
## ----eval=FALSE---------------------------------------------------------------
# load("geno_pheno.RData") # Load data
#
# # Fit model
# fm0 <- fitBLUP(y, K=G)
# fm0$theta <- fm0$varE/fm0$varU # Residual/genetic variances ratio
# fm0$h2 <- fm0$varU/(fm0$varU+fm0$varE) # Heritability
# c(fm0$varU,fm0$varE,fm0$theta,fm0$h2) # Print variance components
#
# save(fm0, file="varComps.RData")
## ----eval=FALSE---------------------------------------------------------------
# nPart <- 5 # Number of partitions
# load("geno_pheno.RData") # Load data
#
# nTST <- ceiling(0.3*length(y)) # Number of elements in TST set
# partitions <- matrix(1,nrow=length(y),ncol=nPart) # Matrix to store partitions
# seeds <- round(seq(1E3, .Machine$integer.max/10, length=nPart))
#
# for(k in 1:nPart){
# set.seed(seeds[k])
# partitions[sample(1:length(y), nTST),k] <- 0
# }
# save(partitions, file="partitions.RData") # Save partitions
## ----eval=FALSE---------------------------------------------------------------
# # Load data
# load("geno_pheno.RData"); load("varComps.RData"); load("partitions.RData")
#
# accSSI <- mu <- varU <- varE <- c() # Objects to store results
#
# for(k in 1:ncol(partitions))
# { cat(" partition = ",k,"\n")
# trn <- which(partitions[,k]==1)
# tst <- which(partitions[,k]==0)
#
# # G-BLUP model
# fm1 <- fitBLUP(y, K=G, trn=trn)
# mu[k] <- fm1$b # Retrieve mu estimate
# varU[k] <- fm1$varU # Retrieve varU
# varE[k] <- fm1$varE # Retrieve varE
#
# # Sparse SI
# fm2 <- SGP(y, K=G, b=mu[k], varU=varU[k], varE=varE[k],
# trn=trn, tst=tst, mc.cores=1, nlambda=100)
# fm3 <- summary(fm2) # Useful function to get results
#
# accuracy <- c(GBLUP=cor(fm1$yHat[tst],y[tst]), fm3$accuracy[,1])/sqrt(fm0$h2)
# lambda <- c(min(fm3$lambda),fm3$lambda[,1])
# nsup <- c(max(fm3$nsup),fm3$nsup[,1])
# accSSI <- rbind(accSSI,data.frame(rep=k,SSI=names(accuracy),accuracy,lambda,nsup))
# }
# save(mu,varU,varE,accSSI,file="results_accuracy.RData")
## ----eval=FALSE---------------------------------------------------------------
# load("results_accuracy.RData")
#
# dat <- data.frame(do.call(rbind,lapply(split(accSSI,accSSI$SSI),
# function(x) apply(x[,-c(1:2)],2,mean))))
# dat$Model <- unlist(lapply(strsplit(rownames(dat),"\\."),function(x)x[1]))
#
# dat2 <- do.call(rbind,lapply(split(dat,dat$Mod),function(x)x[which.max(x$acc),]))
#
# if(requireNamespace("ggplot2")){
# ggplot2::ggplot(dat[dat$nsup>1,],ggplot2::aes(-log(lambda),accuracy)) +
# ggplot2::geom_hline(yintercept=dat["GBLUP",]$accuracy, linetype="dashed") +
# ggplot2::geom_line(ggplot2::aes(color=Model),size=1.1) + ggplot2::theme_bw() +
# ggplot2::geom_point(data=dat2,ggplot2::aes(color=Model),size=2.5)
# }
## ----eval=FALSE---------------------------------------------------------------
# load("geno_pheno.RData"); load("varComps.RData")
# load("partitions.RData"); load("results_accuracy.RData")
#
# lambdaCV <- accSSI_CV <- nsupCV <- c() # Objects to store results
#
# for(k in 1:ncol(partitions))
# { cat(" partition = ",k,"\n")
# trn <- which(partitions[,k]==1)
#
# # Cross-validating the training set
# fm1 <- SGP.CV(y, K=G, trn=trn, nlambda=100, mc.cores=1, nfolds=5, nCV=1)
# lambdaCV[k] <- summary(fm1)$optCOR["lambda"]
#
# # Fit a SSI with the estimated lambda
# fm2 <- SGP(y, K=G, b=mu[k], varU=varU[k], varE=varE[k],
# trn=trn, tst=tst, lambda=lambdaCV[k])
#
# accSSI_CV[k] <- summary(fm2)$accuracy/sqrt(fm0$h2)
# nsupCV <- cbind(nsupCV, fm2$nsup)
# }
# save(accSSI_CV,lambdaCV,nsupCV,file="results_accuracyCV.RData")
## ----eval=FALSE---------------------------------------------------------------
# load("results_accuracy.RData"); load("results_accuracyCV.RData")
#
# dat <- data.frame(GBLUP=accSSI[accSSI$SSI=="GBLUP",]$acc,SSI=accSSI_CV)
# rg <- range(dat)
# tmp <- c(mean(rg),diff(rg)*0.4)
#
# if(requireNamespace("ggplot2")){
# ggplot2::ggplot(dat,ggplot2::aes(GBLUP,SSI)) +
# ggplot2::geom_abline(slope=1,linetype="dotted") + ggplot2::theme_bw() +
# ggplot2::geom_point(shape=21,color="orange") + ggplot2::lims(x=rg,y=rg) +
# ggplot2::annotate("text",tmp[1],tmp[1]-tmp[2],label=round(mean(dat$GBLUP),3)) +
# ggplot2::annotate("text",tmp[1]-tmp[2],tmp[1],label=round(mean(dat$SSI),3))
# }
## ----eval=FALSE---------------------------------------------------------------
# load("results_accuracyCV.RData")
#
# dat <- data.frame(nsup=as.vector(nsupCV))
#
# bw <- round(diff(range(dat$nsup))/40)
# if(requireNamespace("ggplot2")){
# ggplot2::ggplot(data=dat,ggplot2::aes(nsup,stat(count)/length(nsupCV))) +
# ggplot2::theme_bw() +
# ggplot2::geom_histogram(color="gray20",fill="lightblue",binwidth=bw) +
# ggplot2::labs(x=bquote("Support set size(" *n[sup]*")"),y="Frequency")
# }
## ----eval=FALSE---------------------------------------------------------------
# # Load data
# load("geno_pheno.RData"); load("partitions.RData"); load("results_accuracyCV.RData")
#
# part <- 1 # Choose any partition from 1,…,nPart
# trn <- which(partitions[,k]==1)
# tst <- which(partitions[,k]==0)
#
# # Fit SSI with lambda previously estimated using CV
# fm <- SGP(y, K=G, trn=trn, tst=tst, lambda=lambdaCV[part])
#
# plot(net(fm, K=G), i=1:16, unified=FALSE, main=NULL, bg.col="white",
# set.size=c(3,1.5,0.2), point.color="gray40", axis.labels=FALSE)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.