local pls regression for nanoparticles

3 files changed, 33 insertions, 20 deletions

diff --git a/lib/nanoparticle.rb b/lib/nanoparticle.rb
index ca79a3d..65aab23 100644
--- a/lib/nanoparticle.rb
+++ b/lib/nanoparticle.rb
@@ -6,9 +6,10 @@ module OpenTox
     field :core, type: Hash, default: {}
     field :coating, type: Array, default: []
     field :proteomics, type: Hash, default: {}
+
+    attr_accessor :scaled_values
  
     def physchem_neighbors min_sim: 0.9, dataset_id:, prediction_feature_id:
-      p self.name
       dataset = Dataset.find(dataset_id)
       relevant_features = {}
       measurements = []
@@ -52,7 +53,9 @@ module OpenTox
           common_descriptors = relevant_features.keys & substance.physchem_descriptors.keys
           # scale values
           query_descriptors = common_descriptors.collect{|d| (physchem_descriptors[d].median-relevant_features[d]["mean"])/relevant_features[d]["sd"]}
+          @scaled_values = common_descriptors.collect{|d| [d,(physchem_descriptors[d].median-relevant_features[d]["mean"])/relevant_features[d]["sd"]]}.to_h
           neighbor_descriptors = common_descriptors.collect{|d| (substance.physchem_descriptors[d].median-relevant_features[d]["mean"])/relevant_features[d]["sd"]}
+          neighbor_scaled_values = common_descriptors.collect{|d| [d,(substance.physchem_descriptors[d].median-relevant_features[d]["mean"])/relevant_features[d]["sd"]]}.to_h
           #weights = common_descriptors.collect{|d| 1-relevant_features[d]["p_value"]}
           weights = common_descriptors.collect{|d| relevant_features[d]["r"]**2}
           sim = Algorithm::Similarity.weighted_cosine(query_descriptors,neighbor_descriptors,weights)
@@ -61,12 +64,16 @@ module OpenTox
             "measurements" => values,
             "similarity" => sim,
             "common_descriptors" => common_descriptors.collect do |id|
-              {:id => id, :p_value => relevant_features[id]["p_value"], :r_squared => relevant_features[id]["r"]**2}
+              {
+                :id => id,
+                :scaled_value => neighbor_scaled_values[id],
+                :p_value => relevant_features[id]["p_value"],
+                :r_squared => relevant_features[id]["r"]**2}
             end
           } if sim >= min_sim
         end
       end
-      p neighbors.size
+      $logger.debug "#{self.name}: #{neighbors.size} neighbors"
       neighbors.sort!{|a,b| b["similarity"] <=> a["similarity"]}
       neighbors
     end
diff --git a/lib/regression.rb b/lib/regression.rb
index cffcbbf..5028c78 100644
--- a/lib/regression.rb
+++ b/lib/regression.rb
@@ -73,23 +73,19 @@ module OpenTox
 
         activities = []
         weights = []
-        pc_ids = neighbors.collect{|n| Substance.find(n["_id"]).physchem_descriptors.keys}.flatten.uniq
+        pc_ids = neighbors.collect{|n| n["common_descriptors"].collect{|d| d[:id]}}.flatten.uniq.sort
         data_frame = []
         data_frame[0] = []
         
         neighbors.each_with_index do |n,i|
-          neighbor = Substance.find(n["_id"])
-          activities = neighbor["measurements"]
+          activities = n["measurements"]
           activities.each do |act|
             data_frame[0][i] = act
             weights << n["similarity"]
-            neighbor.physchem_descriptors.each do |pid,values| 
-              values = [values] unless values.is_a? Array
-              values.uniq!
-              warn "More than one value for '#{Feature.find(pid).name}': #{values.join(', ')}. Using the median." unless values.size == 1
-              j = pc_ids.index(pid)+1
+            n["common_descriptors"].each do |d| 
+              j = pc_ids.index(d[:id])+1
               data_frame[j] ||= []
-              data_frame[j][i] = values.for_R
+              data_frame[j][i] = d[:scaled_value]
             end
           end if activities
           (0..pc_ids.size+1).each do |j| # for R: fill empty values with NA
@@ -97,10 +93,12 @@ module OpenTox
             data_frame[j][i] ||= "NA"
           end
         end
+
         remove_idx = []
         data_frame.each_with_index do |r,i|
           remove_idx << i if r.uniq.size == 1 # remove properties with a single value
         end
+
         remove_idx.reverse.each do |i|
           data_frame.delete_at i
           pc_ids.delete_at i
@@ -112,7 +110,7 @@ module OpenTox
           prediction
         else
           query_descriptors = pc_ids.collect do |i|
-            substance.physchem_descriptors[i] ? substance.physchem_descriptors[i].for_R : "NA"
+            substance.scaled_values[i] ? substance.scaled_values[i] : "NA"
           end
           remove_idx = []
           query_descriptors.each_with_index do |v,i|
@@ -127,10 +125,9 @@ module OpenTox
           if prediction.nil?
             prediction = local_weighted_average substance, neighbors
             prediction[:warning] = "Could not create local PLS model. Using weighted average of similar substances."
-            prediction
-          else
-            prediction
           end
+          p prediction
+          prediction
         end
       
       end
@@ -172,10 +169,15 @@ rlib = File.expand_path(File.join(File.dirname(__FILE__),"..","R"))
           R.eval "fingerprint <- data.frame(rbind(c(#{query_feature_values.join ','})))"
           R.eval "names(fingerprint) <- features" 
           R.eval "prediction <- predict(model,fingerprint)"
+          value = R.eval("prediction").to_f
+          rmse = R.eval("getTrainPerf(model)$TrainRMSE").to_f
+          r_squared = R.eval("getTrainPerf(model)$TrainRsquared").to_f
+          prediction_interval = value-1.96*rmse, value+1.96*rmse
           {
-            :value => R.eval("prediction").to_f,
-            :rmse => R.eval("getTrainPerf(model)$TrainRMSE").to_f,
-            :r_squared => R.eval("getTrainPerf(model)$TrainRsquared").to_f,
+            :value => value,
+            :rmse => rmse,
+            :r_squared => r_squared,
+            :prediction_interval => prediction_interval
           }
         rescue 
           return nil
diff --git a/lib/validation-statistics.rb b/lib/validation-statistics.rb
index e42d298..6b252b1 100644
--- a/lib/validation-statistics.rb
+++ b/lib/validation-statistics.rb
@@ -100,6 +100,8 @@ module OpenTox
         # TODO: predictions within prediction_interval
         self.rmse = 0
         self.mae = 0
+        #self.within_prediction_interval = 0
+        #self.outside_prediction_interval = 0
         x = []
         y = []
         predictions.each do |cid,pred|
@@ -109,6 +111,9 @@ module OpenTox
             error = pred[:value]-pred[:measurements].median
             self.rmse += error**2
             self.mae += error.abs
+            #if pred[:prediction_interval]
+              #if pred[:measurements]
+            #end
           else
             warnings << "No training activities for #{Compound.find(compound_id).smiles} in training dataset #{model.training_dataset_id}."
             $logger.debug "No training activities for #{Compound.find(compound_id).smiles} in training dataset #{model.training_dataset_id}."
@@ -118,7 +123,6 @@ module OpenTox
         R.assign "prediction", y
         R.eval "r <- cor(measurement,prediction,use='pairwise')"
         self.r_squared = R.eval("r").to_ruby**2
-
         self.mae = self.mae/predictions.size
         self.rmse = Math.sqrt(self.rmse/predictions.size)
         $logger.debug "R^2 #{r_squared}"