demo/fitBLUP.R

In MarcooLopez/SFSI_data: Sparse Family and Selection Index

  setwd(tempdir())
  library(SFSI)
  data(wheatHTP)
  X = scale(X[1:300,])        # Subset and scale markers
  G = tcrossprod(X)/ncol(X)   # Genomic relationship matrix
  y = scale(Y[1:300,"YLD"])   # Subset response variable

  # Fit model, whole data
  fm = fitBLUP(y,K=G)
  fm$varU
  fm$varE
  fm$h2
  cor(y,fm$u)                 # Prediction accuracy

  # Training and testing sets
  tst = sample(seq_along(y),ceiling(0.3*length(y)))
  trn = seq_along(y)[-tst]

  yNA <- y
  yNA[tst] <- NA

  # Fit model, split data
  fm = fitBLUP(yNA,K=G)
  plot(y[tst],fm$u[tst])      # Predicted vs observed values in testing set
  cor(y[tst],fm$u[tst])       # Prediction accuracy in testing set
  cor(y[trn],fm$u[trn])       # Prediction accuracy in training set
  fm$h2                       # Heritability (in training set)