Skip to content

Commit

Permalink
imputeRobust test commented
Browse files Browse the repository at this point in the history
  • Loading branch information
@matthias-da
matthias-da committed Nov 10, 2023
1 parent c6214f0 commit 023300c
Showing 1 changed file with 177 additions and 177 deletions.
354 changes: 177 additions & 177 deletions tests/test_imputeRobust.R
View file
Original file line number Diff line number Diff line change
@@ -1,180 +1,180 @@
library(VIM)
library(laeken)
library(robustbase)
library(mice)
data("Animals_na")

source("~/workspace/VIM/R/imputeRobust.R")

# debugonce(imputeRobust)
imputeRobust(data = Animals_na,
form = formula("lbrain ~ lbody"),
method = "lm", uncert = "wresid", boot = FALSE)

source("~/workspace/vimbook/R/ani_plot.R")

data(Animals, package = "MASS")
ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm", mainlab = NULL,
seed = 1)
ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm-chain", mainlab = NULL,
seed = 1)


ani_plot(Animals = Animals, method = "imputeRobust-lts-normal", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-resid", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-wresid", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm-nboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-normal-nboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-resid-nboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-wresid-nboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm-noboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-normal-noboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-resid-noboot", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lts-wresid-noboot", mainlab = NULL, seed = 123)


ani_plot(Animals = Animals, method = "imputeRobust-gam-pmm", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-gam-normal", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-gam-resid", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-gam-wresid", mainlab = NULL, seed = 123)

ani_plot(Animals = Animals, method = "imputeRobust-lm-pmm", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lm-normal", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lm-resid", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-lm-wresid", mainlab = NULL, seed = 123)

ani_plot(Animals = Animals, method = "imputeRobust-MM-pmm", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-MM-normal", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-MM-resid", mainlab = NULL, seed = 123)
ani_plot(Animals = Animals, method = "imputeRobust-MM-wresid", mainlab = NULL, seed = 123)

debugonce("imputeRobustChain")

Animals_na$lbody[1:3] <- NA
xi <- imputeRobustChain(formulas = list(formula("lbody ~ lbrain"),
formula("lbrain ~ lbody")),
data = Animals_na,
boot = TRUE,
robustboot = TRUE,
method = "lts",
takeAll = TRUE,
eps = 0.5,
maxit = 4,
alpha = 0.5,
uncert = "pmm",
familiy = "Gaussian",
value_back = "matrix",
trace = TRUE)

data(iris)
iris[1:3, 1] <- iris[4:6, 2] <- iris[5:7, 3] <- iris[7:9, 4] <- NA
iris[9:11, 5] <- NA
debugonce("imputeRobustChain")
xi <- imputeRobustChain(formulas = list(formula("Sepal.Length ~ ."),
NULL, NULL, NULL, formula("Species ~ .")),
data = iris,
boot = TRUE,
robustboot = TRUE,
method = "lts",
multinom.method = "multinom",
takeAll = TRUE,
eps = 0.5,
maxit = 4,
alpha = 0.5,
uncert = "pmm",
familiy = "Gaussian",
value_back = "matrix",
trace = TRUE)



## blasting

library(stmo)
data(blasting)
head(blasting)

lm_blasting1 <- lm(log(vibration) ~ log(distance) + log(load),
data = blasting)
summary(lm_blasting1)
lm_blasting2 <- lm(log10(vibration) ~ log10(distance) * location + load,
data = blasting)
summary(lm_blasting2)

blasting$vibration[sample(1:nrow(blasting), 7)] <- NA
w <- complete.cases(blasting)

forms <- list(NULL, NULL, NULL, formula("log10(vibration) ~ log10(distance) * location + load"))
debugonce("imputeRobustChain")
xi <- suppressWarnings(imputeRobustChain(formulas = forms, data = blasting, trace = TRUE, maxit = 1, uncert = "resid"))
xi2 <- suppressWarnings(imputeRobustChain(data = blasting, trace = TRUE, maxit = 20))

plot(xi[, c("vibration", "distance")], col = w + 3)
plot(xi2[, c("vibration", "distance")], col = w + 3)

library(mice)
x2 <- complete(mice(blasting, m = 1))

plot(x2[, c("vibration", "distance")], col = w + 3)


f <- formula("log10(vibration) ~ log10(distance) * location + load")
fn <- f[2]
cn <- "vibration"

data("cherry")
cherry.lm <- lm(log(volume) ~ log(diameter) + log(height), data = cherry)

data("leafburn", package = "faraway")
pairs(leafburn)
# log10(burntimei) = β0 + β1nitrogeni + β2log10(chlorine)i + β3potassiumi +

data(wood, package = "robustbase")
?wood

setMissOut <- function(x, vars = 1:ncol(x), r = rep(0.1, ncol(x))){
# outlier detection
o <- robCompositions::outCoDa(x, coda = FALSE, quantile = 0.99)
n_outliers <- sum(o$outlierIndex)
n_nonoutliers <- sum(!o$outlierIndex)
# sort data
x_na <- x_orig <- rbind(x[!o$outlierIndex, ], x[o$outlierIndex, ])
# set missings in all variables
# but not in outliers
n <- nrow(x)
n_good <- sum(!o$outlierIndex)
n_out <- sum(o$outlierIndex)
for(i in vars){
s <- sample(1:n_good, round(r[i] * n), replace = FALSE)
x_na[s, i] <- NA
}
m1 <- is.na(x_na[, vars])
return(x_na)
}

wood_na <- setMissOut(wood)
xi <- imputeRobustChain(data = wood_na)



setMiss <- function(x, kind = "MCAR", rate = 0.3){
# simple function assuming x has variables Y and X
SEQ <- 1:nrow(x)
n <- nrow(x)
p <- ncol(x)
if(kind == "MCAR"){
s <- sample(SEQ, size = round(rate * n, 0))
}
if(kind == "MAR"){
s <- sample(SEQ, size = round(rate * n, 0), prob = (x$X+abs(min(x$X))+1)/sum(x$X+abs(min(x$X))+1))
}
if(kind == "MNAR"){
s <- sample(SEQ, size = round(rate * n, 0), prob = (x$Y+abs(min(x$Y))+1)/sum(x$Y+abs(min(x$Y))+1))
}
x$Y[s] <- NA
return(x)
}

# library(VIM)
# library(laeken)
# library(robustbase)
# library(mice)
# data("Animals_na")
#
# source("~/workspace/VIM/R/imputeRobust.R")
#
# # debugonce(imputeRobust)
# imputeRobust(data = Animals_na,
# form = formula("lbrain ~ lbody"),
# method = "lm", uncert = "wresid", boot = FALSE)
#
# source("~/workspace/vimbook/R/ani_plot.R")
#
# data(Animals, package = "MASS")
# ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm", mainlab = NULL,
# seed = 1)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm-chain", mainlab = NULL,
# seed = 1)
#
#
# ani_plot(Animals = Animals, method = "imputeRobust-lts-normal", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-resid", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-wresid", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm-nboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-normal-nboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-resid-nboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-wresid-nboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-pmm-noboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-normal-noboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-resid-noboot", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lts-wresid-noboot", mainlab = NULL, seed = 123)
#
#
# ani_plot(Animals = Animals, method = "imputeRobust-gam-pmm", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-gam-normal", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-gam-resid", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-gam-wresid", mainlab = NULL, seed = 123)
#
# ani_plot(Animals = Animals, method = "imputeRobust-lm-pmm", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lm-normal", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lm-resid", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-lm-wresid", mainlab = NULL, seed = 123)
#
# ani_plot(Animals = Animals, method = "imputeRobust-MM-pmm", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-MM-normal", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-MM-resid", mainlab = NULL, seed = 123)
# ani_plot(Animals = Animals, method = "imputeRobust-MM-wresid", mainlab = NULL, seed = 123)
#
# debugonce("imputeRobustChain")
#
# Animals_na$lbody[1:3] <- NA
# xi <- imputeRobustChain(formulas = list(formula("lbody ~ lbrain"),
# formula("lbrain ~ lbody")),
# data = Animals_na,
# boot = TRUE,
# robustboot = TRUE,
# method = "lts",
# takeAll = TRUE,
# eps = 0.5,
# maxit = 4,
# alpha = 0.5,
# uncert = "pmm",
# familiy = "Gaussian",
# value_back = "matrix",
# trace = TRUE)
#
# data(iris)
# iris[1:3, 1] <- iris[4:6, 2] <- iris[5:7, 3] <- iris[7:9, 4] <- NA
# iris[9:11, 5] <- NA
# debugonce("imputeRobustChain")
# xi <- imputeRobustChain(formulas = list(formula("Sepal.Length ~ ."),
# NULL, NULL, NULL, formula("Species ~ .")),
# data = iris,
# boot = TRUE,
# robustboot = TRUE,
# method = "lts",
# multinom.method = "multinom",
# takeAll = TRUE,
# eps = 0.5,
# maxit = 4,
# alpha = 0.5,
# uncert = "pmm",
# familiy = "Gaussian",
# value_back = "matrix",
# trace = TRUE)
#
#
#
# ## blasting
#
# library(stmo)
# data(blasting)
# head(blasting)
#
# lm_blasting1 <- lm(log(vibration) ~ log(distance) + log(load),
# data = blasting)
# summary(lm_blasting1)
# lm_blasting2 <- lm(log10(vibration) ~ log10(distance) * location + load,
# data = blasting)
# summary(lm_blasting2)
#
# blasting$vibration[sample(1:nrow(blasting), 7)] <- NA
# w <- complete.cases(blasting)
#
# forms <- list(NULL, NULL, NULL, formula("log10(vibration) ~ log10(distance) * location + load"))
# debugonce("imputeRobustChain")
# xi <- suppressWarnings(imputeRobustChain(formulas = forms, data = blasting, trace = TRUE, maxit = 1, uncert = "resid"))
# xi2 <- suppressWarnings(imputeRobustChain(data = blasting, trace = TRUE, maxit = 20))
#
# plot(xi[, c("vibration", "distance")], col = w + 3)
# plot(xi2[, c("vibration", "distance")], col = w + 3)
#
# library(mice)
# x2 <- complete(mice(blasting, m = 1))
#
# plot(x2[, c("vibration", "distance")], col = w + 3)
#
#
# f <- formula("log10(vibration) ~ log10(distance) * location + load")
# fn <- f[2]
# cn <- "vibration"
#
# data("cherry")
# cherry.lm <- lm(log(volume) ~ log(diameter) + log(height), data = cherry)
#
# data("leafburn", package = "faraway")
# pairs(leafburn)
# # log10(burntimei) = β0 + β1nitrogeni + β2log10(chlorine)i + β3potassiumi +
#
# data(wood, package = "robustbase")
# ?wood
#
# setMissOut <- function(x, vars = 1:ncol(x), r = rep(0.1, ncol(x))){
# # outlier detection
# o <- robCompositions::outCoDa(x, coda = FALSE, quantile = 0.99)
# n_outliers <- sum(o$outlierIndex)
# n_nonoutliers <- sum(!o$outlierIndex)
# # sort data
# x_na <- x_orig <- rbind(x[!o$outlierIndex, ], x[o$outlierIndex, ])
# # set missings in all variables
# # but not in outliers
# n <- nrow(x)
# n_good <- sum(!o$outlierIndex)
# n_out <- sum(o$outlierIndex)
# for(i in vars){
# s <- sample(1:n_good, round(r[i] * n), replace = FALSE)
# x_na[s, i] <- NA
# }
# m1 <- is.na(x_na[, vars])
# return(x_na)
# }
#
# wood_na <- setMissOut(wood)
# xi <- imputeRobustChain(data = wood_na)
#
#
#
# setMiss <- function(x, kind = "MCAR", rate = 0.3){
# # simple function assuming x has variables Y and X
# SEQ <- 1:nrow(x)
# n <- nrow(x)
# p <- ncol(x)
# if(kind == "MCAR"){
# s <- sample(SEQ, size = round(rate * n, 0))
# }
# if(kind == "MAR"){
# s <- sample(SEQ, size = round(rate * n, 0), prob = (x$X+abs(min(x$X))+1)/sum(x$X+abs(min(x$X))+1))
# }
# if(kind == "MNAR"){
# s <- sample(SEQ, size = round(rate * n, 0), prob = (x$Y+abs(min(x$Y))+1)/sum(x$Y+abs(min(x$Y))+1))
# }
# x$Y[s] <- NA
# return(x)
# }
#

# imputeRobust(x = Animals_na,
# formulas = list(formula("lbody ~ lbrain"),
Expand Down

0 comments on commit 023300c

Please sign in to comment.