# --------------------------------------------- #
# Arabidopsis example #
# --------------------------------------------- #
require(bigRR)
data(Arabidopsis)
X <- matrix(1, length(y), 1)
## Not run:
# # splitting the genotype data into two pieces and re-saving in DatABEL format
# #
# dimnames(Z) <- list(NULL, NULL)
# Z <- scale(Z)
# matrix2databel(Z[,1:100000], 'part1')
# matrix2databel(Z[,100001:ncol(Z)], 'part2')
#
# # fitting SNP-BLUP, i.e. a ridge regression on all the markers across the genome
# #
# SNP.BLUP.result <- hugeRR(y = y, X = X, Z.name = 'part', Z.index = 1:2,
# family = binomial(link = 'logit'), save.cache = TRUE)
#
# # re-run SNP-BLUP - a lot faster since cache data are stored
# SNP.BLUP.result <- hugeRR(y = y, X = X, Z.name = 'part', Z.index = 1:2,
# family = binomial(link = 'logit'))
#
# # fitting HEM, i.e. a generalized ridge regression with marker-specific shrinkage
# #
# HEM.result <- hugeRR_update(SNP.BLUP.result, Z.name = 'part', Z.index = 1:2,
# family = binomial(link = 'logit'))
#
# # plot and compare the estimated effects from both methods
# #
# split.screen(c(1, 2))
# split.screen(c(2, 1), screen = 1)
# screen(3); plot(abs(SNP.BLUP.result$u), cex = .6, col = 'slateblue')
# screen(4); plot(abs(HEM.result$u), cex = .6, col = 'olivedrab')
# screen(2); plot(abs(SNP.BLUP.result$u), abs(HEM.result$u), cex = .6, pch = 19,
# col = 'darkmagenta')
#
# # create a random new genotypes for 10 individuals with the same number of markers
# # and predict the outcome using the fitted HEM
# #
# Z.new <- matrix(sample(c(-1, 1), 10*ncol(Z), TRUE), 10)
# y.predict <- as.numeric(HEM.result$beta + Z.new %*% HEM.result$u)
# #
# # NOTE: The above prediction may not be good due to the scaling in the HEM
# # fitting above, and alternatively, one can either remove the scaling
# # above or scale Z.new by row-binding it with the original Z matrix.
# ## End(Not run)
Run the code above in your browser using DataLab