inst/doc/vigettes_mlm4omics.R

In mlm4omics: Multilevel Model for Multivariate Responses with Missing Values

## ----setup, include=TRUE---------------------------------------------------
library(knitr)
library(rmarkdown)
knitr::opts_chunk$set(echo = TRUE)

## ----installation1,eval=FALSE----------------------------------------------
#  if (!requireNamespace("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("mlm4omics")
#  library(mlm4omics)

## ----installation2,eval=FALSE----------------------------------------------
#  devtools::install_github("superOmics/mlm4omics")
#  library(mlm4omics)

## ----read data,eval=TRUE---------------------------------------------------
library(mlm4omics)
data(pdata,package="mlm4omics")
pdata[1:2,]

## ----table for miss and censor,eval=TRUE-----------------------------------
table(pdata$miss, pdata$censor)

## ----re-assign miss and censor,eval=TRUE-----------------------------------
n=dim(pdata)[1]
for (i in seq_len(n))
if (pdata$miss[i]==1 && pdata$censor[i]==1) pdata$miss[i]=0

## ----set formular examples,eval=TRUE---------------------------------------
formula_completed=var1~var2+treatment;
formula_missing=miss~var2;
formula_censor=censor~1;
formula_subject=~treatment;
response_censorlim=0.002

## ----example 1a,eval=TRUE--------------------------------------------------
model1 <- mlmc(formula_completed = var1~var2+treatment,
formula_missing = miss~var2, 
formula_censor = censor~1, formula_subject = ~sid, 
pdata = pdata, response_censorlim = 0.002, respond_dep_missing = TRUE, 
pidname = "geneid", sidname = "sid", 
iterno = 100, chains = 2, savefile = TRUE)

## ----example 1b-1c,eval=FALSE----------------------------------------------
#  #Example 1b
#  model1 <- mlmc(formula_completed=var1~var2+treatment,
#  formula_missing=miss~var2,
#  formula_censor=censor~1, formula_subject=~sid,
#  pdata=pdata, response_censorlim=0.002,
#  respond_dep_missing=TRUE, pidname="geneid", sidname="sid",
#  alpha_prior <- c(0,0.001), iterno=100, chains=2,
#  savefile=TRUE)
#  
#  #Example 1c
#  prec_example <- matrix(c(0.01,0.001,0.001,0.01),nrow=2,ncol=2)
#  
#  model3 <- mlmc(formula_completed=var1~var2,
#  formula_missing=miss~var2,
#  formula_censor=censor~1, formula_subject=~sid+treatment,
#  pdata=pdata, response_censorlim=0.002, respond_dep_missing=TRUE,
#  pidname="geneid", sidname="sid", prec_prior=prec_example,
#  iterno=1000, chains=2, savefile=TRUE)
#  

## ----example 2a,eval=FALSE-------------------------------------------------
#  model5 <- mlmm(formula_completed=var1~var2+treatment,
#  formula_missing=miss~var2,
#  formula_subject=~sid, pdata=pdata,
#  respond_dep_missing=FALSE, pidname="geneid", sidname="sid",
#  iterno=1000, chains=2, savefile=FALSE)

## ----example 2b prior,eval=TRUE--------------------------------------------
prec_example <- matrix(c(0.01,0.001,0.001,0.001,0.01,0.001,0.001,0.001,0.01),
nrow=3, ncol=3)

## ----example 2b, 2c,eval=FALSE---------------------------------------------
#  #Example 2b use both priors
#  model4 <- mlmm(formula_completed=var1~var2+treatment,
#  formula_missing=miss~var2,
#  formula_subject=~sid, pdata=pdata,
#  respond_dep_missing=FALSE, pidname="geneid", sidname="sid",
#  alpha_prior <- c(0,0.001), prec_prior=prec_example,
#  iterno=100, chains=2, savefile=TRUE)
#  
#  #Example 2c. Use alpha prior
#  model5 <- mlmm(formula_completed=var1~var2+treatment,
#  formula_missing=miss~var2,
#  formula_subject=~sid+treatment, pdata=pdata,
#  respond_dep_missing=FALSE, pidname="geneid",alpha_prior <- c(0,0.001),
#  sidname="sid", iterno=100, chains=2, savefile=TRUE)
#  

## ----plot one,echo=TRUE,eval=TRUE------------------------------------------
summaryreader <- read.csv(file = file.path(getwd(),"outsummary.csv"),
header=TRUE, sep=",", skip=0)

iterno <- dim(summaryreader)[1]; burnin <- iterno/2

U.1.1 <- rowMeans(matrix(c(summaryreader$chain.1.U.1.1,
summaryreader$chain.2.U.1.1), nrow=iterno, ncol=2))[burnin:iterno]

meanU <- mean(U.1.1)
qU <- quantile(U.1.1,p <- seq(0, 1, by=0.025))
scale <- seq(0, 1, by=0.025)

plot(scale, qU, pch=19, ylab="quantiles of estimate", xlab="quantiles")

segments(0,qU[names(qU)=="50%"],1,qU[names(qU)=="50%"],lwd=2,col="red")
segments(0,qU[names(qU)=="2.5%"],1,qU[names(qU)=="2.5%"],lty=2,lwd=2,col="red")
segments(0,qU[names(qU)=="97.5%"],1,qU[names(qU)=="97.5%"],lty=2,lwd=2,
col="red")

legend(0.5,qU[names(qU)=="50%"],"median",cex=0.8,bty="n")
legend(0.03,qU[names(qU)=="2.5%"],"2.5%",cex=0.8,bty="n")
legend(0.90,qU[names(qU)=="97.5%"],"97.5%",cex=0.8,bty="n")
qU

## ----plot two,echo=TRUE,eval=TRUE------------------------------------------
sample1reader <- read.csv(file <- file.path(getwd(),"samples_1.csv"),
header=TRUE, sep=",", skip=25)

sample2reader <- read.csv(file <- file.path(getwd(),"samples_2.csv"),
header=TRUE, sep=",", skip=25)

#plot variable U.1.1 - the intercept of first unit
trajectory_length <- dim(sample1reader)[1]

plot(seq(1, trajectory_length, by=1), sample1reader$U.1.1, xlab="trajectory
number", ylab="U.1.1", type="n",
ylim = c(min(sample1reader$U.1.1, sample2reader$U.1.1, na.rm=TRUE),
max(sample1reader$U.1.1, sample2reader$U.1.1, na.rm=TRUE)))

trajectory <- seq(1, trajectory_length, by=1)

lines(trajectory, sample1reader$U.1.1)
lines(trajectory, sample2reader$U.1.1, col="red")

## ----sessionInfo,echo=TRUE,eval=TRUE---------------------------------------
sessionInfo()