Revision a53065a09c3fce65a63e137deb5bccb6162e6cff authored by Matthias Templ on 18 November 2020, 20:10:02 UTC, committed by cran-robot on 18 November 2020, 20:10:02 UTC
1 parent 9ae1e67
Raw File
rz.R
# require("robCompositions")
# require("MASS")
# require("robustbase")
# genVarsSimple <- function(n=200, p=50, k=3){
#   T <- matrix(rnorm(n*k,0,1), ncol=k)
#   #  B <- matrix(rnorm(p*k,0,0.1), ncol=k)
#   B <- matrix(runif(p*k,-1,1), ncol=k)
#   X <- T %*% t(B)
#   E <-  matrix(rnorm(n*p, 0,0.1), ncol=p)
#   XE <- X + E
#   XE <- isomLRinv(XE)
#   return(constSum(XE))  
# }
# 
# # 
# # crnorm <- function(n, mu, Sigma){ 
# #   constSum(data.frame(exp(mvrnorm(n, mu, Sigma))))
# # }
# # sigGen <- function(p, d){
# #   x <- diag(p)
# #   x[upper.tri(x)] <- x[lower.tri(x)] <- d
# #   x
# # }
# set.seed(1234)
# # x <- crnorm(50, rep(10,10), Sigma=sigGen(10,0.9))
# # x <- xOrig <- x+rnorm(ncol(x)*nrow(x), 10, 1)
# x <- xOrig <- data.frame(genVarsSimple(50, k=3, p=12))
# 
# xOrig <- x
# lim <- apply(x, 2, quantile, 0.05)
# for(i in 1:ncol(x)){
#   x[x[,i] < lim[i],i] <- 0
# }
# xOrig[1,1] <- xOrig[2,2] <- xOrig[3,1] <- 0.00001
# x[1,1] <- x[2,2] <- x[3,1] <- 0
# 
# 
# res01 <- imputeBDLs(x, dl=dl, method="pls", variation = FALSE, verbose = TRUE)
# res02 <- imputeBDLs(x, dl=dl, method="lm", variation = TRUE, verbose = TRUE)
# res01b <- imputeBDLs(x, dl=dl, method="pls", variation = FALSE, verbose = TRUE, test = TRUE)
# res02b <- imputeBDLs(x, dl=dl, method="lm", variation = TRUE, verbose = TRUE, test = TRUE)
# 
# 
# res1 <- impRZilr(x, dl=dl, method="lm")
# res2 <- impRZilr(x, dl=dl, method="MM")
# res3 <- impRZilr(x, dl=dl, method="pls", nComp="boot", verbose=TRUE)
# res4 <- impRZilr(x, dl=dl, method="pls", nComp=rep(5,ncol(x)), verbose=TRUE)
# 
# epsilon <- 55*.Machine$double.eps
# if(any(!c(
#   isTRUE(is.logical((x[1,3]/x[1,4] - res1$x[1,3]/res1$x[1,4]) < epsilon)),
#   isTRUE(is.logical((x[1,3]/x[1,4] - res2$x[1,3]/res2$x[1,4]) < epsilon)),
#   isTRUE(is.logical((x[1,3]/x[1,4] - res3$x[1,3]/res3$x[1,4]) < epsilon)),
#   isTRUE(is.logical((x[1,3]/x[1,4] - res4$x[1,3]/res4$x[1,4]) < epsilon)),
#   isTRUE(is.logical((x[1,3]/x[1,4] - res01$x[1,3]/res01$x[1,4]) < epsilon)),
#   isTRUE(is.logical((x[1,3]/x[1,4] - res02$x[1,3]/res01$x[1,4]) < epsilon))
# ))){
#   stop("ratios are not preserved.")
# }
# 
# check1 <- function(x, w, dl){
#   check <- logical(ncol(x))
#   for(i in 1:ncol(x)){
#     check[i] <- any(x[w[,i],i] > dl[i])
#   }
#   if(any(check)){
#     res <- FALSE
#     stop("values are over detection limit")
#   } else {
#     res <- TRUE
#     message("imputed data passed check below detection limit")
#   }
#   invisible(res)
# }
# 
# if(all(!c(
#   check1(res1$x, w, dl=dl),
#   check1(res2$x, w, dl=dl),
#   check1(res3$x, w, dl=dl),
#   check1(res4$x, w, dl=dl),
#   check1(res01$x, w, dl=dl),
#   check1(res02$x, w, dl=dl)
#   ))){
#     stop("one method imputed above detection limit")
# }
# 
# ### quality
# require("matrixcalc")
# rdcm <- function(x, y){
#   ocov <- cov(isomLR(x))
#   rcov <- cov(isomLR(y))
#   return(frobenius.norm(ocov-rcov)/frobenius.norm(ocov))
# }
# 
# ced <- function(x, y, ni){
#   return(aDist(x, y)/ni)
# }
# 
# ni <- sum(x == 0)
# 
# ced(res01$x, xOrig, ni = ni)
# ced(res02$x, xOrig, ni = ni)
# ced(res01b$x, xOrig, ni = ni)
# ced(res02b$x, xOrig, ni = ni)
# ced(res1$x, xOrig, ni = ni)
# ced(res2$x, xOrig, ni = ni)
# ced(res3$x, xOrig, ni = ni)
# ced(res4$x, xOrig, ni = ni)
# rdcm(res01$x, xOrig)
# rdcm(res02$x, xOrig)
# rdcm(res01b$x, xOrig)
# rdcm(res02b$x, xOrig)
# rdcm(res1$x, xOrig)
# rdcm(res2$x, xOrig)
# rdcm(res3$x, xOrig)
# rdcm(res4$x, xOrig)
# 
# ## compare with zCompositions:
# require("zCompositions")
# resz1 <- multRepl(x, label = 0, dl = dl)
# zm1 <- ced(resz1, xOrig, ni = ni)
# zm1
# zm2 <- rdcm(resz1, xOrig)
# zm2
# zm1 > ced(res01$x, xOrig, ni = ni)
# zm1 > ced(res02$x, xOrig, ni = ni)
# zm2 > rdcm(res01$x, xOrig)
# zm2 > rdcm(res02$x, xOrig)
# 
# 
# ## HD data:
# set.seed(1234)
# x <- xOrig <- crnorm(50, rep(10,60), Sigma=sigGen(60,0.9))
# lim <- 0.01
# x[x < lim] <- 0
# w <- x==0
# dl <- rep(lim, ncol(x))
# 
# resHD01 <- imputeBDLs(x, dl=dl, method="pls", nComp="boot", verbose=TRUE, variation = FALSE, maxit=5)
# resHD02 <- imputeBDLs(x, dl=dl, method="pls", nComp=rep(5,ncol(x)), verbose=TRUE, , maxit=5, variation = FALSE)
# resHD01var <- imputeBDLs(x, dl=dl, method="lm", nComp="boot", verbose=TRUE, variation = TRUE, maxit=5)
# resHD01b <- imputeBDLs(x, dl=dl, method="pls", nComp="boot", verbose=TRUE, variation = FALSE, test=TRUE, maxit=5)
# resHD02b <- imputeBDLs(x, dl=dl, method="pls", nComp=rep(5,ncol(x)), verbose=TRUE, maxit=5, variation = FALSE, test=TRUE)
# resHD01bvar <- imputeBDLs(x, dl=dl, method="lm", nComp="boot", verbose=TRUE, variation = TRUE, test=TRUE, maxit=5)
# 
# 
# resHD1 <- impRZilr(x, dl=dl, method="pls", nComp="boot", verbose=TRUE, maxit=5)
# resHD2 <- impRZilr(x, dl=dl, method="pls", nComp=rep(5,ncol(x)), verbose=TRUE, maxit=5)
# 
# ni <- sum(x == 0)
# 
# ced(resHD01$x, xOrig, ni = ni)
# ced(resHD02$x, xOrig, ni = ni)
# ced(resHD01var$x, xOrig, ni = ni)
# ced(resHD01b$x, xOrig, ni = ni)
# ced(resHD02b$x, xOrig, ni = ni)
# ced(resHD01bvar$x, xOrig, ni = ni)
# ced(resHD1$x, xOrig, ni = ni)
# ced(resHD2$x, xOrig, ni = ni)
# 
# rdcm(resHD01$x, xOrig)
# rdcm(resHD02$x, xOrig)
# rdcm(resHD01var$x, xOrig)
# rdcm(resHD01b$x, xOrig)
# rdcm(resHD02b$x, xOrig)
# rdcm(resHD01bvar$x, xOrig)
# rdcm(resHD1$x, xOrig)
# rdcm(resHD2$x, xOrig)
# 
# ## compare with zCompositions:
# require("zCompositions")
# resz1 <- multRepl(x, label = 0, dl = dl)
# zm1 <- ced(resz1, xOrig, ni = ni)
# zm1
# zm2 <- rdcm(resz1, xOrig)
# zm2
# zm1 > ced(resHD01$x, xOrig, ni = ni)
# zm1 > ced(resHD02$x, xOrig, ni = ni)
# zm2 > rdcm(resHD01$x, xOrig)
# zm2 > rdcm(resHD02$x, xOrig)
# 
# 
# # 
# # 
# # data(arcticLake)
# # x <- arcticLake
# # ## generate rounded zeros artificially:
# # x[x[,1] < 10, 1] <- 0
# # xia <- impRZilr(x, dl=c(10,44,0), eps=0.01, method="MM")
# # xia$x
# 
# 
# #data(expenditures)
# #xOrig <- x <- expenditures
# 
# ## DL as negative values
# #x[x < 350] <- - 350
# 
# #impCoda2(x)
# 
# ## DL given:
# #impCoda2(xOrig, dl=c(0, 350, 350, 350, 350))
# 
# 
# #imp <- impCoda(x, method='roundedZero')
# #imp2 <- alrEM(x, pos=2, dl=rep(5,3))$xImp
# 
# 
back to top