rstudio/keras3#778 https://stackoverflow.com/questions/57242141/how-to-share-weight-between-two-keras-layers https://towardsdatascience.com/build-the-right-autoencoder-tune-and-optimize-using-pca-principles-part-ii-24b9cca69bd6
visualizating with tensorboard https://medium.com/@xianbao.qian/keras-model-visualization-on-tensorboard-8da47a6442bf
Create a model and write graph to logdir folder import tensorflow as tf graph = tf.Graph() m = build_model() # Your model implementation with graph.as_default():
m.compile(loss=’categorical_crossentropy’, optimizer=’adam’, metrics=[‘accuracy’]) writer = tf.summary.FileWriter(logdir=’logdir’, graph=graph) writer.flush() Then start TensorBoard on that folder from console tensorboard –logdir logdir
Run tensorboard –logdir=YOUR_LOG_DIR –host=127.0.0.1 in command prompt, and type localhost:6006 in chrome, this works for me (Win10, anaconda4.3.16, python3.5.3, tensorflow1.1.0).
Compile the model model.compile(loss=keras.losses.categorical_crossentropy, optimizer=keras.optimizers.Adam(), metrics=[‘accuracy’])
tensorboard = TensorBoard( log_dir=’.\logs’, histogram_freq=1, write_images=True ) keras_callbacks = [ tensorboard ]
model.fit(input_train, target_train, batch_size=batch_size, epochs=no_epochs, verbose=verbosity, validation_split=validation_split, callbacks=keras_callbacks)
source("LatentConfounderBNlearnv2.R")
reticulate::conda_list()
reticulate::use_condaenv("tensorflowv1")
reticulate::source_python("aupy.py")
setwd("/Users/fred/Documents/projects/latvar/")
load("/Users/fred/Documents/projects/latvar/res_alldata.RData", verbose = T)
load("final_model_nolvp_novp.RData", verbose = T)
load('res_missing.RData', verbose = T)
## library(keras)
## keras::is_keras_available()
tensorflow::tf_version()
K <- backend()
K$clear_session()
train = datalist$data_noisy
trainlv = train %>% dplyr::select(-U1.out, -U2.out)
vars = getDrivers(res_alldata, 'Z', maxpath = 2, cutoff = 0.1)$Drivers
vars = getDrivers(res_missing, 'Z', maxpath = 2, cutoff = 0.1)$Drivers
resList = getGraphResiduals (res_missing,
vars,
trainlv)
resDev = residualDeviance(trainlv,
resList, isOrdinal = T)
tst = findLatentVars(resDev,
method = 'autoencoder'
)
plot_history(tst$details$pcObj$history)
tensorboard = TensorBoard(
log_dir='./logs',
histogram_freq=1,
write_images=TRUE
)
K = import("keras")
ModelCheckpoint = K$callbacks$ModelCheckpoint
keras_callbacks = list(
tensorboard
)
mc = ModelCheckpoint("bestmodel.h5", monitor = "val_mean_squared_error", mode = "max",
verbose=1, save_best_only=TRUE)
keras_callbacks = list(
mc
)
tst$details$pcObj$save("tstmodel.h5")
tst = findLatentVars(resDev,
method = 'autoencoder',
callback = keras_callbacks
)
model_visualize = import("plotnn.visualize.model_visualize")
saved_model = load_model('bestmodel.h5')
#Build your model here
model_visualize(model, name="graph")
tf = import("tensorflow")
graph = tf$Graph()
writer = tf$summary$FileWriter("logs", graph)
writer$flush()
K = import("keras")
plot_model = K$utils$vis_utils$plot_model
plot_model(tst$details$pcObj, to_file = "model_ae2.pdf", show_shapes = TRUE, show_layer_names = TRUE)
annv = import("ann_visualizer.visualize")
ann_viz = annv$ann_viz
ann_viz(tst$details$pcObj, filename = 'plotae.gv')
tst$confounders %>%
cor(train$U2.out)
tst$confounders %>%
cor(train$U1.out)
testlin = findLatentVars(resDev,
method = 'linear'
)
testlin$confounders %>%
cor(train$U2.out)
testlin$confounders %>%
cor(train$U1.out)
myplot = function(tstlin, train){
toplot = cbind(as.data.frame(tstlin$confounder), train[, c("U1.out", "U2.out")]) %>% gather(var, val, -U1.out, -U2.out)
ggp1 = toplot%>% ggplot(aes(x = U1.out, y = val)) + geom_point() + facet_wrap( ~ var)
ggp2 = toplot %>% ggplot(aes(x = U2.out, y = val)) + geom_point() + facet_wrap( ~ var)
cowplot::plot_grid(ggp1, ggp2, ncol = 1)
}
##myplot(tst, train)
myplot(tstlin, train)
tstlin$confounders %>%
cor(train$U1.out)
tstlin$confounders %>%
cor(train$U2.out)
(arch = rev(round(exp(seq(log(4), log(min(200, ncol(resDev)/2)), length.out=4)))))
fwrite(resDev, file = "resDev.csv")
tstaelin = findLatentVars(resDev,
scale. = T,
method = 'autoencoder',
architecture=2,
activation = 'linear',
drRate=0.0,
use_batch_norm=F,
nIter = 500,
batch_size = 32,
optimizer = "RMSprop",
metrics = 'mse', learning_rate = 0.05,
fname = "fitae_1layer_pca_rmsprop_lr005_itera500.pdf"
)
tstaelin = findLatentVars(resDev,
scale. = T,
method = 'autoencoder',
architecture=c(20, 10, 5, 3),
activation = 'relu',
drRate=0.3,
use_batch_norm=F,
nIter = 500,
batch_size = 32,
optimizer = "RMSprop",
metrics = 'mse', learning_rate = 0.05,
fname = "fitae_rmsprop_lr005_itera500.pdf"
)
tstaelin = findLatentVars(resDev,
scale. = T,
method = 'autoencoder',
architecture=NULL,
activation = 'sigmoid',
activation_coding = "sigmoid",
activation_output = "linear",
drRate=0.2,
use_batch_norm=T,
nIter =100,
batch_size = 16,
optimizer = "RMSprop",
metrics = 'mse',
fname = "fitae_rmsprop_lr005_itera500_new.pdf"
)
plot_history(tstaelin$details$pcObj$history)
qplot(tstaelin$confounders[[1]])
qplot(tstaelin$confounders[[2]])
tstaelin$confounders %>%
cor(train$U1.out)
tstaelin$confounders %>%
cor(train$U2.out)
myplot(tstaelin, train)
newaepred = fread("predae.csv", data.table = F)
newaepred %>% cor(train$U1.out)
newaepred %>% cor(train$U2.out)
setwd("/Users/fred/Documents/projects/latvar/")
latvarae = readRDS("aetest/latVars.RDS")source("LatentConfounderBNlearn.R")
load("final_model_nolvp_novp.RData", verbose = T)
train = datalist$data_noisy
trainlv = train %>% dplyr::select(-U1.out, -U2.out)
blacklistlv = rbind(data.frame(from = "Z", to = colnames(trainlv)))
library(doParallel)
cl <- makeCluster(5) ## for multi-threading
registerDoParallel(cl)
res_missing_small = getEnsemble2(trainlv, blacklist = blacklistlv,
restart = 100, Nboot = 10,
prior = "vsp",
score = "bge",
algorithm = 'hc',
parallel = TRUE
)
test_wrong = latentDiscovery(
res_missing_small,
nItera=5,
data = trainlv,
"Z",
workpath="pca_wrong",
freqCutoff = 0.01,
maxpath = 1,
alpha = 0.01,
scale. = TRUE,
method = "linear",
latent_iterations = 100,
truecoef = datalist$coef %>% filter(output=="Z"),
truelatent=train %>% dplyr::select("U1.out","U2.out"),
include_downstream = TRUE,
include_output = TRUE,
multiple_comparison_correction = T,
debug = F,
parallel = TRUE
)
source("LatentConfounderBNlearnv2.R")
test_right = latentDiscovery(
res_missing_small,
nItera=5,
data = trainlv,
"Z",
workpath="pca_right",
freqCutoff = 0.01,
maxpath = 1,
alpha = 0.01,
scale. = TRUE,
method = "linear",
latent_iterations = 100,
truecoef = datalist$coef %>% filter(output=="Z"),
truelatent=train %>% dplyr::select("U1.out","U2.out"),
include_downstream = TRUE,
include_output = TRUE,
multiple_comparison_correction = T,
debug = F,
parallel = TRUE,
wrongway = TRUE ## this undo the fix in getGraphResiduals
)library(keras)
DenseTied <- R6::R6Class("DenseTied",
inherit = KerasLayer,
public = list(
master_layer = NULL,
output_dim = NULL,
weights = NULL,
bias = NULL,
initialize = function(master_layer = NULL) {
self$master_layer = master_layer
},
build = function(input_shape) {
self$weights = k_transpose(self$master_layer$weights[[1]])
self$output_dim <- self$weights$shape$as_list()[[2]]
self$bias <- self$add_weight(
name = 'bias',
shape = list(self$output_dim),
initializer = initializer_constant(0),
trainable = TRUE
)
message("build worked fine")
},
call = function(x, mask = NULL) {
message("in call")
browser()
res = k_dot(x, self$weights) + self$bias
message("finished call")
res
},
compute_output_shape = function(input_shape) {
message("in shape")
res = list(input_shape[[1]], self$output_dim)
message("finished shape")
res
}
)
)
layer_densetied <- function(object, master_layer, name = NULL, trainable = TRUE) {
browser()
create_layer(DenseTied, object, list(
master_layer = master_layer,
name = name,
trainable = trainable
)
)
}
input_layer <- layer_input(shape = ncol(mtcars))
l1 = layer_dense(units = 100, input_shape = 11)
l2 = layer_dense(units = 10)
l3 = layer_dense(units = 2)
decl1 = layer_tied_dense(master_layer = l3)
decl2 = layer_tied_dense(master_layer = l2)
decl3 = layer_tied_dense(master_layer = l1)
output = input_layer %>%
l1() %>%
layer_activation('relu') %>%
layer_dropout(0.2) %>%
l2() %>%
layer_activation("relu") %>%
layer_dropout(0.2) %>%
l3() %>%
layer_activation("relu") %>%
decl1() %>%
layer_activation("relu") %>%
decl2() %>%
layer_activation("relu") %>%
decl3() %>%
layer_activation("linear")
model <- keras_model(inputs = input_layer,
outputs = output
)
model %>%
compile(
loss = "mse",
optimizer = "adam"
)
history = keras::fit(model, as.matrix(mtcars),
as.matrix(mtcars),
epochs=1000,
batch_size=32,
shuffle=TRUE,
##validation_data= list(x_train, x_train)
validation_split = 0.1
)
plot(history)
model$layers[[5]]$weights[[1]] %>% dim
model$layers[[10]]$get_weights()
input <- layer_input(shape = ncol(mtcars))
dense_1 <- layer_dense(units = 128)
dense_2 <- layer_dense(units = 256)
dense_3 <- layer_dense(units = 256)
dense_6 <- layer_dense(units = 128)
dense_5 <- layer_dense(units = 256)
dense_4 <- layer_dense(units = 256)
out = layer_dense(units = 11)
dense_1_transposed <- layer_tied_dense(master_layer = dense_1)
dense_2_transposed <- layer_tied_dense(master_layer = dense_2)
dense_3_transposed <- layer_tied_dense(master_layer = dense_3)
output <- input %>%
dense_1() %>%
layer_activation("selu") %>%
dense_2() %>%
layer_activation("selu") %>%
dense_3() %>%
layer_activation("selu") %>%
layer_dropout(0.2) %>%
dense_3_transposed() %>%
layer_activation("selu") %>%
dense_2_transposed() %>%
layer_activation("selu") %>%
dense_1_transposed() %>%
layer_activation("selu")
model %>%
compile(
loss = "mse",
optimizer = "adam"
)
keras::fit(model, as.matrix(mtcars),
as.matrix(mtcars),
epochs=10,
batch_size=16,
shuffle=TRUE,
##validation_data= list(x_train, x_train)
validation_split = 0.1
)
output <- input %>%
dense_1() %>%
layer_activation("relu") %>%
dense_1_transposed() %>%
layer_activation("relu")
input <- layer_input(shape = ncol(mtcars))
output <- input %>%
dense_1() %>%
layer_activation("selu") %>%
dense_2() %>%
layer_activation("selu") %>%
dense_3() %>%
layer_activation("selu") %>%
layer_dropout(0.65) %>%
dense_4() %>%
layer_activation("selu") %>%
dense_5() %>%
layer_activation("selu") %>%
dense_6() %>%
layer_activation("selu") %>%
out()
masked_mse <- function(y_true, y_pred) {
mask_true <- k_cast(k_not_equal(y_true, 0), k_floatx())
masked_squared_error <- k_square(mask_true * (y_true - y_pred))
masked_mse <- k_sum(masked_squared_error)/k_sum(mask_true)
masked_mse
}
rmse <- function(y_true, y_pred) {
masked_mse(y_true, y_pred) ^ 0.5
}
model %>%
compile(
loss = "mse",
metrics = list(rmse = rmse),
optimizer = "adam"
)
keras::fit(model, as.matrix(mtcars),
as.matrix(mtcars),
epochs=10,
batch_size=16,
shuffle=TRUE,
##validation_data= list(x_train, x_train)
validation_split = 0.1
)
model %>%
fit_generator(
sparse_generator(as.matrix(mtcars), 128),
epochs = 100,
steps_per_epoch = nrow(as.matrix(mtcars))/128,
callbacks = callback_tensorboard()
)
evaluate_generator(model, sparse_generator(netflix3m$test, batch_size = 128), steps = 1000)
split_ind <- iris$Species %>% caret::createDataPartition(p = 0.8,list = FALSE)
train <- iris[split_ind,]
test <- iris[-split_ind,]
train_X <- train[,1:4] %>% as.matrix()
train_y <- train[,5] %>%as.integer %>%
keras::to_categorical()
test_X <- test[,1:4] %>% as.matrix()
input_layer <-
layer_input(shape = c(4))
encoder <-
input_layer %>%
layer_dense(units = 150, activation = "relu") %>%
layer_batch_normalization() %>%
layer_dropout(rate = 0.2) %>%
layer_dense(units = 50, activation = "relu") %>%
layer_dropout(rate = 0.1) %>%
layer_dense(units = 25, activation = "relu") %>%
layer_dense(units = 2) # 2 dimensions for the output layer
decoder <-
encoder %>%
layer_dense(units = 150, activation = "relu") %>%
layer_dropout(rate = 0.2) %>%
layer_dense(units = 50, activation = "relu") %>%
layer_dropout(rate = 0.1) %>%
layer_dense(units = 25, activation = "relu") %>%
layer_dense(units = 4) # 4 dimensions for the original 4 variables
autoencoder_model <- keras_model(inputs = input_layer, outputs = decoder)
autoencoder_model %>% compile(
loss='mean_squared_error',
optimizer='adam',
metrics = c('accuracy')
)
summary(autoencoder_model)
history <-
autoencoder_model %>%
keras::fit(train_X,
train_X,
epochs=100,
shuffle=TRUE,
validation_data= list(test_X, test_X)
)
input_layer <-
layer_input(shape = c(4))
dense_1 <- layer_dense(units = 150, activation = 'relu')
dense_2 <- layer_dense(units = 3, activation = "relu")
dense_2_t = layer_tied_dense(master_layer = dense_2)
dense_1_t = layer_tied_dense(master_layer = dense_1)
encoder <-
input_layer %>%
dense_1() %>%
dense_2_t()
decoder <-
encoder %>%
dense_2_t() %>%
dense_1_t()
autoencoder_model <- keras_model(inputs = input_layer, outputs = decoder)
autoencoder_model %>% compile(
loss='mean_squared_error',
optimizer='adam',
metrics = c('accuracy')
)
summary(autoencoder_model)
history <-
autoencoder_model %>%
keras::fit(train_X,
train_X,
epochs=100,
shuffle=TRUE,
validation_data= list(test_X, test_X)
)
## try this one
dense_1 <- layer_dense(units = 128)
dense_2 <- layer_dense(units = 10)
dense_3 <- layer_dense(units = 2)
dense_1_transposed <- layer_tied_dense(master_layer = dense_1)
dense_2_transposed <- layer_tied_dense(master_layer = dense_2)
dense_3_transposed <- layer_tied_dense(master_layer = dense_3)
input <- layer_input(shape = ncol(mtcars))
output <- input %>%
dense_1() %>%
layer_activation("selu") %>%
dense_2() %>%
layer_activation("selu") %>%
dense_3() %>%
layer_activation("selu") %>%
layer_dropout(0.2) %>%
dense_3_transposed() %>%
layer_activation("selu") %>%
dense_2_transposed() %>%
layer_activation("selu") %>%
dense_1_transposed() %>%
layer_activation("selu")
model <- keras_model(input, output)
model %>%
compile(
loss = "mse",
optimizer = "adam"
)
keras::fit(model, as.matrix(mtcars),
as.matrix(mtcars),
epochs=10,
batch_size=16,
shuffle=TRUE,
##validation_data= list(x_train, x_train)
validation_split = 0.1
)
model$layers[[4]]$kernel
model$layers[[9]]$kernel
model$layers[[9]]$get_weights()
dense_1 <- layer_dense(units = 128)
dense_2 <- layer_dense(units = 10)
dense_3 <- layer_dense(units = 2)
dense_1_transposed <- DenseTiedLayer$new(tied_to = dense_1, units = ncol(128))
dense_2_transposed <- DenseTiedLayer$new(tied_to = dense_2, units = 10)
dense_3_transposed <- DenseTiedLayer$new(tied_to = dense_3, units = 2)
input <- layer_input(shape = ncol(mtcars))
output <- input %>%
dense_1() %>%
layer_activation("selu") %>%
dense_2() %>%
layer_activation("selu") %>%
dense_3() %>%
layer_activation("selu") %>%
layer_dropout(0.2) %>%
dense_3_transposed() %>%
layer_activation("selu") %>%
dense_2_transposed() %>%
layer_activation("selu") %>%
dense_1_transposed() %>%
layer_activation("selu")
model <- keras_model(input, output)
model %>%
compile(
loss = "mse",
optimizer = "adam"
)
keras::fit(model, as.matrix(mtcars),
as.matrix(mtcars),
epochs=10,
batch_size=16,
shuffle=TRUE,
##validation_data= list(x_train, x_train)
validation_split = 0.1
)
model$layers[[4]]$kernel
model$layers[[9]]$kernelhttps://stackoverflow.com/questions/53751024/tying-autoencoder-weights-in-a-dense-keras-layer
https://stackoverflow.com/questions/53751024/tying-autoencoder-weights-in-a-dense-keras-layer
class DenseTied(Layer):
def __init__(self, units,
activation=None,
use_bias=True,
kernel_initializer='glorot_uniform',
bias_initializer='zeros',
kernel_regularizer=None,
bias_regularizer=None,
activity_regularizer=None,
kernel_constraint=None,
bias_constraint=None,
tied_to=None,
**kwargs):
self.tied_to = tied_to
if 'input_shape' not in kwargs and 'input_dim' in kwargs:
kwargs['input_shape'] = (kwargs.pop('input_dim'),)
super().__init__(**kwargs)
self.units = units
self.activation = activations.get(activation)
self.use_bias = use_bias
self.kernel_initializer = initializers.get(kernel_initializer)
self.bias_initializer = initializers.get(bias_initializer)
self.kernel_regularizer = regularizers.get(kernel_regularizer)
self.bias_regularizer = regularizers.get(bias_regularizer)
self.activity_regularizer = regularizers.get(activity_regularizer)
self.kernel_constraint = constraints.get(kernel_constraint)
self.bias_constraint = constraints.get(bias_constraint)
self.input_spec = InputSpec(min_ndim=2)
self.supports_masking = True
def build(self, input_shape):
assert len(input_shape) >= 2
input_dim = input_shape[-1]
if self.tied_to is not None:
self.kernel = K.transpose(self.tied_to.kernel)
self._non_trainable_weights.append(self.kernel)
else:
self.kernel = self.add_weight(shape=(input_dim, self.units),
initializer=self.kernel_initializer,
name='kernel',
regularizer=self.kernel_regularizer,
constraint=self.kernel_constraint)
if self.use_bias:
self.bias = self.add_weight(shape=(self.units,),
initializer=self.bias_initializer,
name='bias',
regularizer=self.bias_regularizer,
constraint=self.bias_constraint)
else:
self.bias = None
self.input_spec = InputSpec(min_ndim=2, axes={-1: input_dim})
self.built = True
def compute_output_shape(self, input_shape):
assert input_shape and len(input_shape) >= 2
output_shape = list(input_shape)
output_shape[-1] = self.units
return tuple(output_shape)
def call(self, inputs):
output = K.dot(inputs, self.kernel)
if self.use_bias:
output = K.bias_add(output, self.bias, data_format='channels_last')
if self.activation is not None:
output = self.activation(output)
return outputfrom numpy.random import seed
seed(123)
from tensorflow import set_random_seed
##from tensorflow.compat.v1 import set_random_seed
set_random_seed(234)
import sklearn
from sklearn import datasets
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn import decomposition
import scipy
import tensorflow as tf
from keras.models import Model, load_model
from keras.layers import Input, Dense, Layer, InputSpec
from keras.callbacks import ModelCheckpoint, TensorBoard
from keras import regularizers, activations, initializers, constraints, Sequential
from keras import backend as K
from keras.constraints import UnitNorm, Constraint
import pandas as pd
exec(open("aupy.py").read())
df = pd.read_csv("/Users/fred/Documents/projects/latvar/resDev.csv", header = 0)
X_train, X_test = train_test_split(df, test_size=0.5, random_state=123)
# Scale the data between 0 and 1.
scaler = MinMaxScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
X_test_scaled = scaler.transform(X_test)
X_train_scaled
scaler2=MinMaxScaler()
scaler2.fit(df)
df_scaled=scaler2.transform(df)
nb_epoch = 500
batch_size = 16
input_dim = X_train_scaled.shape[1] #num of predictor variables,
encoding_dim = 2
learning_rate = 1e-3
encoder = Dense(encoding_dim, activation="linear", input_shape=(input_dim,), use_bias = True, kernel_regularizer=WeightsOrthogonalityConstraint(encoding_dim, weightage=1., axis=0), kernel_constraint=UnitNorm(axis=0))
decoder = DenseTied(input_dim, activation="linear", tied_to=encoder, use_bias = False)
autoencoder = Sequential()
autoencoder.add(encoder)
autoencoder.add(decoder)
autoencoder.compile(metrics=['accuracy'],
loss='mean_squared_error',
optimizer='sgd')
autoencoder.summary()
autoencoder.fit(X_train_scaled, X_train_scaled,
epochs=nb_epoch,
batch_size=batch_size,
shuffle=True,
verbose=0)
train_predictions = autoencoder.predict(X_train_scaled)
print('Train reconstrunction error\n', sklearn.metrics.mean_squared_error(X_train_scaled, train_predictions))
test_predictions = autoencoder.predict(X_test_scaled)
print('Test reconstrunction error\n', sklearn.metrics.mean_squared_error(X_test_scaled, test_predictions))
latvar=encoder.predict(X_train_scaled)
train_predictions = autoencoder.predict(X_train_scaled)
autoencoder.fit(df_scaled,df_scaled,
epochs=nb_epoch,
batch_size=batch_size,
shuffle=True,
verbose=0)
enc=Sequential()
enc.add(encoder)
from keras.utils import plot_model
plot_model(autoencoder, to_file='model.png')
train_predictions_enc = enc.predict(df_scaled)
np.savetxt("predae.csv", train_predictions_enc, delimiter=",")https://stackoverflow.com/questions/53751024/tying-autoencoder-weights-in-a-dense-keras-layer
https://stackoverflow.com/questions/53751024/tying-autoencoder-weights-in-a-dense-keras-layer
from numpy.random import seed
seed(123)
import sklearn
from sklearn import datasets
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn import decomposition
import scipy
import tensorflow as tf
from keras.models import Model, load_model
from keras.layers import Input, Dense, Layer, InputSpec
from keras.callbacks import ModelCheckpoint, TensorBoard
from keras import regularizers, activations, initializers, constraints, Sequential
from keras import backend as K
from keras.constraints import UnitNorm, Constraint
import pandas as pd
class DenseTied(Layer):
def __init__(self, units,
activation=None,
use_bias=True,
kernel_initializer='glorot_uniform',
bias_initializer='zeros',
kernel_regularizer=None,
bias_regularizer=None,
activity_regularizer=None,
kernel_constraint=None,
bias_constraint=None,
tied_to=None,
**kwargs):
self.tied_to = tied_to
if 'input_shape' not in kwargs and 'input_dim' in kwargs:
kwargs['input_shape'] = (kwargs.pop('input_dim'),)
super().__init__(**kwargs)
self.units = units
self.activation = activations.get(activation)
self.use_bias = use_bias
self.kernel_initializer = initializers.get(kernel_initializer)
self.bias_initializer = initializers.get(bias_initializer)
self.kernel_regularizer = regularizers.get(kernel_regularizer)
self.bias_regularizer = regularizers.get(bias_regularizer)
self.activity_regularizer = regularizers.get(activity_regularizer)
self.kernel_constraint = constraints.get(kernel_constraint)
self.bias_constraint = constraints.get(bias_constraint)
self.input_spec = InputSpec(min_ndim=2)
self.supports_masking = True
def build(self, input_shape):
assert len(input_shape) >= 2
input_dim = input_shape[-1]
if self.tied_to is not None:
self.kernel = K.transpose(self.tied_to.kernel)
self._non_trainable_weights.append(self.kernel)
else:
self.kernel = self.add_weight(shape=(input_dim, self.units),
initializer=self.kernel_initializer,
name='kernel',
regularizer=self.kernel_regularizer,
constraint=self.kernel_constraint)
if self.use_bias:
self.bias = self.add_weight(shape=(self.units,),
initializer=self.bias_initializer,
name='bias',
regularizer=self.bias_regularizer,
constraint=self.bias_constraint)
else:
self.bias = None
self.input_spec = InputSpec(min_ndim=2, axes={-1: input_dim})
self.built = True
def compute_output_shape(self, input_shape):
assert input_shape and len(input_shape) >= 2
output_shape = list(input_shape)
output_shape[-1] = self.units
return tuple(output_shape)
def call(self, inputs):
output = K.dot(inputs, self.kernel)
if self.use_bias:
output = K.bias_add(output, self.bias, data_format='channels_last')
if self.activation is not None:
output = self.activation(output)
return output
class WeightsOrthogonalityConstraint (Constraint):
def __init__(self, encoding_dim, weightage = 1.0, axis = 0):
self.encoding_dim = encoding_dim
self.weightage = weightage
self.axis = axis
def weights_orthogonality(self, w):
if(self.axis==1):
w = K.transpose(w)
if(self.encoding_dim > 1):
m = K.dot(K.transpose(w), w) - K.eye(self.encoding_dim)
return self.weightage * K.sqrt(K.sum(K.square(m)))
else:
m = K.sum(w ** 2) - 1.
return m
def __call__(self, w):
return self.weights_orthogonality(w)
class UncorrelatedFeaturesConstraint (Constraint):
def __init__(self, encoding_dim, weightage = 1.0):
self.encoding_dim = encoding_dim
self.weightage = weightage
def get_covariance(self, x):
x_centered_list = []
for i in range(self.encoding_dim):
x_centered_list.append(x[:, i] - K.mean(x[:, i]))
x_centered = tf.stack(x_centered_list)
covariance = K.dot(x_centered, K.transpose(x_centered)) / tf.cast(x_centered.get_shape()[0], tf.float32)
return covariance
# Constraint penalty
def uncorrelated_feature(self, x):
if(self.encoding_dim <= 1):
return 0.0
else:
output = K.sum(K.square(
self.covariance - tf.math.multiply(self.covariance, K.eye(self.encoding_dim))))
return output
def __call__(self, x):
self.covariance = self.get_covariance(x)
return self.weightage * self.uncorrelated_feature(x)
library(reticulate)
library(tidyverse)
source_python("aupy.py")
sk = import("sklearn", convert = FALSE)
keras = import("keras", convert = FALSE)
K = keras$backend
library(data.table)
df = fread("/Users/fred/Documents/projects/latvar/resDev.csv", data.table = F)
minmax = sk$preprocessing$MinMaxScaler
scaler = minmax()
scaler$fit(df)
df_scaled = scaler$transform(df)
nb_epoch = 500L
batch_size = 16L
input_dim = df_scaled$shape[1] #num of predictor variables,
encoding_dim = 2L
learning_rate = 1e-3
shap0 = 10L
shap1 = 5L
encoder0 = Dense(shap0,
activation = 'relu',
input_shape = list(input_dim),
use_bias = TRUE
)
encoder1 = Dense(shap1,
activation = 'relu',
input_shape = list(shap0),
use_bias = TRUE
)
encoder = Dense(encoding_dim,
activation = 'linear',
input_shape = list(shap1),
use_bias = TRUE,
kernel_regularizer=WeightsOrthogonalityConstraint(encoding_dim, weightage=1., axis=0L),
kernel_constraint=UnitNorm(axis=0L),
activity_regularizer=UncorrelatedFeaturesConstraint(encoding_dim, weightage = 1.)
)
decoder = DenseTied(shap1, activation="relu", tied_to=encoder, use_bias = TRUE)
decoder1 = DenseTied(shap0, activation="relu", tied_to=encoder1, use_bias = TRUE)
decoder0 = DenseTied(input_dim, activation="linear", tied_to=encoder0, use_bias = TRUE)
autoencoder = Sequential()
autoencoder$add(encoder0)
autoencoder$add(encoder1)
autoencoder$add(encoder)
autoencoder$add(decoder)
autoencoder$add(decoder1)
autoencoder$add(decoder0)
autoencoder$compile(metrics=list('mse'),
loss="mean_squared_error",
optimizer='adam')
autoencoder$summary()
history = autoencoder$fit(df_scaled, df_scaled,
epochs=as.integer(1000),
batch_size=16L,
validation_split = 0.1,
shuffle=TRUE
)
plot_history(history,"mse")
res %>% names
res$vas_loss
lossdf = tibble(val_loss = as.numeric(res$val_loss),
loss = as.numeric(res$loss),
Epoch = 1:length(loss)
)
gather(lossdf, key, value, -Epoch) %>% ggplot(aes(x = Epoch, y = value, colour = key)) + geom_point() + geom_line()
enc=Sequential()
enc$add(encoder0)
enc$add(encoder1)
enc$add(encoder)
latevar = enc$predict(df_scaled)
cor(latevar, train$U1.out)
cor(latevar, train$U2.out)
qplot(latevar[, 1])
qplot(latevar[, 2])
qplot(latevar[, 1], train$U2.out)+ ggtitle("cor=", signif(cor(latevar[, 1], train$U1.out), 2))
qplot(latevar[, 1], train$U2.out)+ ggtitle("cor=", signif(cor(latevar[, 1], train$U2.out), 2))
qplot(latevar[, 2], train$U2.out)+ ggtitle("cor=", signif(cor(latevar[, 2], train$U2.out), 2))
qplot(latevar[, 2], train$U1.out) + ggtitle("cor=", signif(cor(latevar[, 2], train$U1.out), 2))
DenseTiedLayer <- R6::R6Class(
"DenseTiedLayer",
inherit = KerasLayer,
public = list(
tied_to = NULL,
units = NULL,
activation = NULL,
use_bias = NULL,
kernel_initializer = NULL,
bias_initializer = NULL,
activity_regularizer = NULL,
kernel_constraint = NULL,
bias_constraint = NULL,
initialize = function(units,
activation = NULL,
use_bias = TRUE,
kernel_initializer = "glorot_uniform",
bias_initializer = "zeros",
bias_regularizer = NULL,
activity_regularizer = NULL,
kernel_constraint = NULL,
bias_constraint = NULL,
tied_to = NULL,
...
) {
self$tied_to <- tied_to
self$units = units
self$activation = activation
self$use_bias = use_bias
self$kernel_initializer = kernel_initializer
self$bias_initializer = bias_initializer
self$activity_regularizer = activity_regularizer
self$kernel_constraint = kernel_constraint
self$bias_constraint = bias_constraint
},
build = function(input_shape) {
input_dim = input_shape[2]
if(!is.null(self$tied_to)){
self$kernel = k_transpose(self$tied_to$kernel)
}else{
self$kernel <- self$add_weight(
name = 'kernel',
shape = list(self$output_dim),
initializer = self$kernel_initializer,
regularizer = self$kernel_regularizer,
constraint = self$kernel_constraint,
trainable = TRUE
)
}
if(self$use_bias){
self$bias <- self$add_weight(
name = 'bias',
shape = list(self$output_dim),
initializer = initializer_constant(0),
regularizer = self$bias_regularizer,
constraint = self$bias_constraint,
trainable = TRUE
)
}else
self$bias = NULL
},
compute_output_shape = function(input_shape) {
list(input_shape[[1]], self$output_dim)
},
call = function(x, mask = NULL) {
output = k_dot(x, self$kernel)
if(self$use_bias)
output = k_bias_add(output, self$bias)
if(!is.null(self$activation)){
output = self$activation(output)
}
return(output)
}
)
)
library(keras)
CustomLayer <- R6::R6Class("CustomLayer",
inherit = KerasLayer,
public = list(
output_dim = NULL,
kernel = NULL,
initialize = function(output_dim) {
self$output_dim <- output_dim
},
build = function(input_shape) {
self$kernel <- self$add_weight(
name = 'kernel',
shape = list(input_shape[[2]], self$output_dim),
initializer = initializer_random_normal(),
trainable = TRUE
)
},
call = function(x, mask = NULL) {
k_dot(x, self$kernel)
},
compute_output_shape = function(input_shape) {
list(input_shape[[1]], self$output_dim)
}
)
)
layer_custom <- function(object, output_dim, name = NULL, trainable = TRUE) {
create_layer(CustomLayer, object, list(
output_dim = as.integer(output_dim),
name = name,
trainable = trainable
))
}
# use it in a model
model <- keras_model_sequential()
model %>%
layer_dense(units = 32, input_shape = c(32,32)) %>%
layer_custom(output_dim = 32)use metric that can be applied to partially directed graph.
- also we may want to update the algorithm to use the sampling in bnlearn to generate the predictions rather than the current way we are using.
THe ideas was to use bnlearn::CPDAG to convert estimated bnlearn object and then to score it. We can then compare it to the output of FCI.
to score the functions instead of a structural score. What is CPDAG?
CPDAG simply compare a dag into a cpdag. this useful to compare resulits of FCI with a dag Some references: https://www.researchgate.net/post/How-can-I-measure-similarity-between-two-networks
shd calculates de hamming distance after converting the networks fo PDAG. This can be used to compare the output of FCI with a dag. THe problem is that we have additional nodes representing the estimated latent variables so we can’t use this function as is.
Maybe the way around this is to calculate the error with respect to the true by removing all edges involving the latent variable. We can compare FCI and the Score based method with respect to the true. If the latent variable is well estimated then you epxcect the skeleton will be closer to the true.
library(bnlearn)
dag1 = model2network("[A][B|A][C|A]")
dag2 = model2network("[A|B:C][B][C]")
dag3 = model2network("[A][B|A]")
dag4 = model2network("[A][B|A][C|A]")
compare(dag1, dag2)
compare(dag3, dag4, arcs = TRUE)
hamming(dag1, dag2)
hamming(dag3, dag4)
cpdag(dag1)
cpdag(dag2)
shd(dag1, dag2)
shd(dag3, dag4)
par(mfrow = c(3, 2))
graphviz.compare(dag1, dag2, dag3, dag4-)Lets a synthetic dataset and estimate the inferred network with a score-based method like tabu search. Assume that the resulting network is true.
library(bnlearn)
library(tidyverse)
data(gaussian.test)
truenet = tabu(gaussian.test) ## true network
res_rm = remove.node(res, "A")
plot(res_rm)
## suppose A is the latent variable
datalat = gaussian.test %>% select(-A)
ddim = dim(datalat)
noisesd =4
datalatnoise = datalat + matrix(rnorm(prod(ddim), sd = noisesd), nrow = ddim[1])
reslat = tabu(datalatnoise)
plot(reslat)
reslat2 = fast.iamb(datalatnoise)
plot(reslat2)
pcalg::fci(datalatnoise)
bnlearn::shd(reslat, res_rm)
bnlearn::shd(reslat2, res_rm)
library(CompareCausalNetworks)
rfcinetmat = getParents(datalatnoise, method = "rfci")
colnames(rfcinetmat) = colnames(datalatnoise)
rownames(rfcinetmat) = colnames(datalatnoise)
rfci_ig = igraph::graph_from_adjacency_matrix(rfcinetmat)
rfci_bn = as.bn(rfci_ig)
plot(rfci_bn)
bnlearn::shd(rfci_bn, truenet_rm)plot(truenet)
Now suppose the node A is latent. The best structure that we can generate would be what we would get by removing node A from the true network:
truenet_rm = remove.node(truenet, "A")
plot(truenet_rm)
Now we can try to learn networks from the dataset with the latent using both score and constraint-based methods.
We remove the node A from the dataset and add noise
datalat = gaussian.test %>% select(-A)
ddim = dim(datalat)
noisesd =4
datalatnoise = datalat + matrix(rnorm(prod(ddim), sd = noisesd), nrow = ddim[1])
Using Tabu:
estnet_tabu = tabu(datalatnoise)
plot(estnet_tabu)
Score for tabu search:
bnlearn::shd(estnet_tabu, truenet_rm)Using Fast Incremental Association (Fast-IAMB) Learning Algorithm:
estnet_iamb = fast.iamb(datalatnoise)
plot(estnet_iamb)
bnlearn::shd(estnet_iamb, truenet_rm)Finally using RFCI
library(CompareCausalNetworks)
rfcinetmat = getParents(datalatnoise, method = "rfci")
colnames(rfcinetmat) = colnames(datalatnoise)
rownames(rfcinetmat) = colnames(datalatnoise)
rfci_ig = igraph::graph_from_adjacency_matrix(rfcinetmat)
rfci_bn = as.bn(rfci_ig)
plot(rfci_bn)
bnlearn::shd(rfci_bn, truenet_rm)one with R code so we can easily test it. Also one that it is easy to use.
Look at some of the papers and how they test their algorithms. See if there is something we can use easily.
- It would be good to have an example showing the advantage of using the residuals over not using them.
devtools::load_all("NetRes")
library(bnlearn)
network = SFNetwork$new(numVertices=50, topology='star')
generated = network$generateData(numSamples = 100, latIdx = 1)
print(network)
plot(network)
bnlearn::graphviz.plot(generated$graph, layout='dot', highlight = list(nodes = grep('U_', colnames(generated$data), value=T), col = "tomato", fill = "orange"))
tabuArgs = list(start = NULL, whitelist = NULL, blacklist = NULL, score = 'bge', prior = 'vsp', debug = FALSE, tabu = 10, max.iter = Inf, maxp = Inf, optimized = TRUE)
##tabuArgs = list(start = NULL, whitelist = NULL, blacklist = NULL, score = 'bge', debug = FALSE, tabu = 10, max.iter = Inf, maxp = Inf, optimized = TRUE)
netResObj = NetRes$new(generated$data, true.graph = generated$graph, nIter = 2, nBoot=100, algorithm='tabu',
algorithm.args = tabuArgs, mode='foo')
##netResObj$assess() #true network has 0 edges which causes crashes
##after inference
bnlearn::graphviz.plot(netResObj$ensemble[[1]][[1]], layout='dot')
bnlearn::graphviz.plot(netResObj$ensemble[[2]][[1]], layout='dot', highlight = list(nodes = grep('U_', bnlearn::nodes(netResObj$ensemble[[2]][[1]]), value=T), col = "tomato", fill = "orange"))
netResObj_oracle = NetRes$new(generated$data, true.graph = generated$graph, nIter = 2, nBoot=100, algorithm='tabu',
algorithm.args = tabuArgs, mode="oracular")
netResObj_oracle2 = NetRes$new(generated$data %>%
mutate(vlat=U_v50) %>%
select(-U_v50),
true.graph = generated$graph, nIter = 2, nBoot=100, algorithm='tabu',
algorithm.args = tabuArgs)
netResObj_oracle2 %>% namessetwd("/Users/fred/Documents/projects/netres_version2")
devtools::load_all("NetRes")
library(bnlearn)
network = SFNetwork$new(numVertices=200)
generated = network$generateData(numSamples=100)
print(network)
plot(network)
bnlearn::graphviz.plot(generated$graph, layout='dot', highlight = list(nodes = grep('U_', colnames(generated$data), value=T), col = "tomato", fill = "orange"))
plotIgraph(as.igraph(generated$graph),nodesep=0.0001,fill=list("U_.*"="orange"))
## check pca
pcaval=PCAtools::pca(generated$data)
corrplot(cbind(pcaval$loadings[,1:2],select(generated$data,matches("U_.*"))))
library(ggplot2)
qplot(pcaval$loadings[,1],generated$data$U_v17)+geom_smooth()
qplot(pcaval$loadings[,2],generated$data$U_v17)+geom_smooth()
qplot(pcaval$loadings[,1],generated$data$U_v80)+geom_smooth()
qplot(pcaval$loadings[,2],generated$data$U_v80)+geom_smooth()
generated$graph$nodes %>% names %>%grep(pattern="U_.*",value=T)
(latvars = colnames(generated$data) %>% grep(pattern="U_.*",value=T))
tabuArgs = list(start = NULL, whitelist = NULL, blacklist = NULL, score = 'bge', prior = 'vsp', debug = FALSE, tabu = 100, max.iter = Inf, maxp = Inf, optimized = TRUE)
netResObj = NetRes$new(generated$data, true.graph = generated$graph, nIter = 3, nBoot=100, algorithm='tabu', algorithm.args = tabuArgs, mode=NULL)
## learn with true latent in there. no latent variables
tstdata=generated$data
ii=which(colnames(tstdata)%in%latvars)
for(jj in 1:length(ii)){
tstdata[[paste0("vlat",jj)]]=tstdata[,ii[jj]]
}
netResObj_nolat = NetRes$new(tstdata,
true.graph = generated$graph, nIter = 3, nBoot=100, algorithm='tabu', algorithm.args = tabuArgs, mode=NULL)
dev.new()
plotIgraph(as.igraph(netResObj_nolat$ensemble[[1]][[2]]),nodesep=0.001,fill=list("vlat.*"='orange',".*PC.*"='yellow'))
## with blacklist
varsnou=grep("^v.*",colnames(generated$data),value=T)
blist = rbind(data.frame(from = varsnou, to = "vlat1"),
data.frame(from = varsnou, to = "vlat2")
)
tabuArgs2 = list(start = NULL, whitelist = NULL, blacklist = blist, score = 'bge', prior = 'vsp', debug = FALSE, tabu = 100, max.iter = Inf, maxp = Inf, optimized = TRUE)
netResObj_nolat = NetRes$new(tstdata,
true.graph = generated$graph, nIter = 3, nBoot=100, algorithm='tabu', algorithm.args = tabuArgs2, mode=NULL,debug=F)
plotIgraph(as.igraph(netResObj_nolat$ensemble[[3]][[3]]),nodesep=0.001,fill=list("vlat.*"='orange',".*PC.*"='yellow'))##library(conflicted)
library(tidyverse)
library(igraph)
library(bnlearn)
library(simcausal)
setwd("/Users/fred/Documents/projects/netres_fork/")
##source("~/gitRepos/gnsutils/gnsutils/R/networks.R")
source("~/Documents/projects/gnsutils/gnsutils/R/networks.R")
##library(brms)
rskew_normal=brms::rskew_normal
devtools::load_all("./NetRes")
##insurancenet=readRDS("~/gitRepos/general-functions/org/bnlearn/insurance.rds")
insurancenet=readRDS("~/Documents/gitRepos/general-functions/org/bnlearn/insurance.rds")
##conflicts_prefer(bnlearn::as.igraph)
insurance_ig=as.igraph(insurancenet)
ppPlotIgraph(insurance_ig,nodesep=0.01)
insurance_ig=addConfounderIgraph(insurance_ig,prob=0.8,ucoef=4)
ppPlotIgraph(insurance_ig,nodesep=0.01)
insurance_simc=igraph2simcausal(insurance_ig,specific_parameters = attributes(insurance_ig)$ucoef)
insurance_simc_set=set.DAG(insurance_simc)
data_insurance=sim(insurance_simc_set,n=1000)
true.graph=bnlearn::as.bn(insurance_ig)
algArgs = list(start = NULL, whitelist = NULL, blacklist = NULL, score = 'ebic-g', prior = 'vsp', debug = FALSE, tabu = 100, max.iter = Inf, maxp = Inf, optimized = TRUE, max.rank=10)
## run with all variables
netResObj_oracle = NetRes$new(
data_insurance %>%
select(-ID),
true.graph = true.graph,
nIter = 1,
nBoot=30,
algorithm='tabu',
algorithm.args = algArgs,
mode='oracular',
weightedResiduals=FALSE,
scale=TRUE,
latentSpaceMethod='pca',
lvPrefix="U",
nCores=4,
#BPPARAM=BiocParallel::SerialParam(),
debug=T)
dev.new()
netResObj_oracle$assess(iteration=1,lvPrefix="U",ci=T)
ppPlotIgraph(as.igraph(netResObj_oracle$ensemble[[1]][[1]]),nodesep=0.01)
## run removing U
netResObj_naive = NetRes$new(
data_insurance %>%
select(-ID),
true.graph = true.graph,
nIter = 1,
nBoot=30,
algorithm='tabu',
algorithm.args = algArgs,
mode='normal',
weightedResiduals=FALSE,
scale=TRUE,
latentSpaceMethod='pca',
lvPrefix="U",
nCores=4,
#BPPARAM=BiocParallel::SerialParam(),
debug=F)
ppPlotIgraph(as.igraph(netResObj_naive$ensemble[[1]][[1]]),nodesep=0.01)
ppPlotIgraph(netResObj_oracle$true.graph %>% as.igraph(),nodesep=0.01)
netResObj_naive$assess(iteration=1,lvPrefix="U",ci=T)
## run with estimated LV
netResObj = NetRes$new(
data_insurance %>%
select(-ID),
true.graph = true.graph,
nIter = 10,
nBoot=100,
algorithm='tabu',
algorithm.args = algArgs,
mode='normal',
weightedResiduals=FALSE,
scale=TRUE,
latentSpaceMethod='pca',
lvPrefix="U",
nCores=4,
#BPPARAM=BiocParallel::SerialParam(),
debug=F)
netResObj$assess(iteration="all")
titer=5
dev.new();netResObj$assess(iteration=titer)
netResObj$assess(iteration=titer,oracle=netResObj_oracle)
netResObj$assess(nCores=4,lvPrefix="U",ci=T,save_to_pdf="insurance.pdf",iteration="all",oracle=netResObj_oracle)
nCores = detectCores()-8
cluster = makeCluster(nCores)
registerDoParallel(cluster)
dud = function(x) x
ens = bn.boot(data_insurance %>%
select(-ID), statistic = dud, R=30, algorithm = "tabu", algorithm.args = algArgs, cluster=cluster)
stopCluster(cluster)
edges=bnlearn::custom.strength(ens,colnames(data_insurance %>%
select(-ID)))
##colnames(true.graph$arcs) =c("from","to")
pred = bnlearn::as.prediction(edges, true.graph)
perf = ROCR::performance(pred, "tpr", "fpr")
(auc = round(ROCR::performance(pred, "auc")@y.values[[1]], 2))
plot(perf, main = paste("AUC:", auc), colorize=TRUE)
(aucpr = round(ROCR::performance(pred, "aucpr")@y.values[[1]], 2))
## with minet
true_adj=igraph::get.adjacency(insurance_ig) %>% as.matrix
est_adj=edges %>% mutate(pred=strength*direction) %>% select(from,to,pred) %>%
igraph::graph_from_data_frame() %>%
igraph::get.adjacency(attr="pred") %>%
as.matrix()
est_adj=est_adj[colnames(true_adj),colnames(true_adj)]
library(minet)
compa=validate(est_adj,true_adj)
auc.roc(compa)
auc.pr(compa)
## with latent variable
cluster = makeCluster(nCores)
registerDoParallel(cluster)
dud = function(x) x
ens_lat = bn.boot(data_insurance %>%
select(-ID,-U), statistic = dud, R=30, algorithm = "tabu", algorithm.args = algArgs, cluster=cluster)
stopCluster(cluster)
edges_lat=bnlearn::custom.strength(ens_lat,colnames(data_insurance %>%
select(-ID,-U)))
## with minet
true_adj_removeu=igraph::get.adjacency(delete_vertices(insurance_ig,"U")) %>% as.matrix
est_adj_lat=edges_lat %>% mutate(pred=strength*direction) %>% select(from,to,pred) %>%
igraph::graph_from_data_frame() %>%
igraph::get.adjacency(attr="pred") %>%
as.matrix()
est_adj_lat=est_adj_lat[colnames(true_adj_removeu),colnames(true_adj_removeu)]
library(minet)
compa_lat=validate(est_adj_lat,true_adj_removeu)
auc.roc(compa_lat)
auc.pr(compa_lat)
##
netResObj = NetRes$new(
data_insurance %>%
select(-ID),
true.graph = true.graph,
nIter = 3,
nBoot=10,
algorithm='tabu',
algorithm.args = algArgs,
mode='normal',
weightedResiduals=FALSE,
scale=TRUE,
latentSpaceMethod='pca',
lvPrefix="U",
nCores=4,
#BPPARAM=BiocParallel::SerialParam(),
debug=F)
netResObj$assess(nCores=4,lvPrefix="U",iteration=2,ci=T)
netResObj$assess(nCores=4,lvPrefix="U",save_to_pdf="testassess.pdf",ci=T)
netResObj_oracle = NetRes$new(
data_insurance %>%
select(-ID),
true.graph = true.graph,
nIter = 2,
nBoot=30,
algorithm='tabu',
algorithm.args = algArgs,
mode='oracular',
weightedResiduals=FALSE,
scale=TRUE,
latentSpaceMethod='pca',
lvPrefix="U",
nCores=4,
#BPPARAM=BiocParallel::SerialParam(),
debug=F)
netResObj_oracle$assess(iteration=1,lvPrefix="U",ci=T)
netResObj$assess(nCores=4,lvPrefix="U",oracle=netResObj_oracle,ci=T,iteration=3)
netResObj$assess(nCores=4,lvPrefix="U",ci=T,iteration='all',save_to_pdf="testall.pdf")set.seed(42)
devtools::load_all("NetRes")
##devtools::load_all("../netres_boris/NetRes/")
library(bnlearn)
network_star = SFNetwork$new(numVertices=50, topology='star', addOneEdge = FALSE)
generated_star = network_star$generateData(numSamples = 100, latIdx = 1)
bnlearn::graphviz.plot(generated_star$graph, layout='dot', highlight = list(nodes = grep('U_', colnames(generated_star$data), value=T), col = "tomato", fill = "orange"))
tabuArgs = list(start = NULL, whitelist = NULL, blacklist = NULL, score = 'ebic-g', debug = FALSE, tabu = 10, max.iter = Inf, maxp = Inf, optimized = TRUE)
library(corrplot)
library(ggplot2)
netResObj_star = NetRes$new(generated_star$data,
true.graph = generated_star$graph, nIter = 5, nBoot=100, algorithm='tabu',
algorithm.args = tabuArgs, mode='foo')
netResObj_star$assess(iteration=2)
## assessment with oracle should not crash but it will just ignore it. Not point of wating effort on this since this is a edge case.
netResObj_star_oracle = NetRes$new(generated_star$data,
true.graph = generated_star$graph, nIter = 5, nBoot=100, algorithm='tabu',
algorithm.args = tabuArgs, mode='oracular')
netResObj_star$assess(iteration=1,lvPrefix='U_.*',oracle=netResObj_star_oracle)- How does the performance compare to simply calculating the PCA on the data and taking the first PC as latent variable and forget about iterations
- what about using the PA to determine the number of PC to use
- what is the effect of getting the wrong number of PC in the estimate
- what if there are no latent confounders. how does the algorithm behave?
notion?
lift curves.
have an option where it plots the result for a specific interation