In stephenslab/flashr: Empirical Bayes Matrix Factorization

knitr::opts_chunk$set(echo = TRUE)

A simple rank 1 simulation

First we will simulate some rank 1 data:

  library("flashr")
  set.seed(1)
  n=100
  p=1000
  ll = rnorm(n)
  ff = rnorm(p)
  LF = outer(ll,ff)
  Y = LF + rnorm(n*p)

Now add some missing data at 50 percent of entries

  Y.miss = Y
  for(i in 1:n){ # set half of Y to be missign at random
    Y.miss[i,sample(1:p,p/2)]=NA
  }

Run flash (with $K=1$) on the full data

  data = flash_set_data(Y)
  f= flash(data,Kmax=1,verbose=TRUE)
  plot(LF,flash_get_fitted_values(f),main="True LF vs fitted value (full data)",xlab="true LF",ylab="fitted LF")
  abline(a=0,b=1,col="red")

And on missing data

  data.miss = flash_set_data(Y.miss)
  f.miss = flash(data.miss,Kmax=1,verbose=TRUE)
  plot(LF,flash_get_fitted_values(f.miss),main="True LF vs fitted values (missing data)",xlab="true LF",ylab="fitted LF")
  abline(a=0,b=1,col="red")

And compute overall RMSE of estimated low-rank structure.

  sqrt(mean((LF-flash_get_fitted_values(f))^2))
  sqrt(mean((LF-flash_get_fitted_values(f.miss))^2))

And try setting the var_type slightly differently:

  f= flash(data,Kmax=1,var_type = "constant",verbose=TRUE)
  plot(LF,flash_get_fitted_values(f),main="True LF vs fitted value (full data)",xlab="true LF",ylab="fitted LF")
  abline(a=0,b=1,col="red")
  sqrt(mean((LF-flash_get_fitted_values(f))^2))

  f.miss = flash(data.miss,Kmax=1,var_type = "constant",verbose=TRUE)
  plot(LF,flash_get_fitted_values(f.miss),main="True LF vs fitted values (missing data)",xlab="true LF",ylab="fitted LF")
  abline(a=0,b=1,col="red")
  sqrt(mean((LF-flash_get_fitted_values(f.miss))^2))

stephenslab/flashr documentation built on Jan. 31, 2024, 2:07 a.m.