[Mattice-commits] r69 - pkg/R

[noreply at r-forge.r-project.org]

Author: andrew_hipp
Date: 2008-12-02 23:31:58 +0100 (Tue, 02 Dec 2008)
New Revision: 69

Modified:
   pkg/R/informationCriterion.R
   pkg/R/summarizingAnalyses.R
Log:
summary now correctly calculates the node probability tree by tree, then averages over all trees as well as conditioned on trees that posses the branch

Modified: pkg/R/informationCriterion.R
===================================================================
--- pkg/R/informationCriterion.R	2008-12-02 21:43:17 UTC (rev 68)
+++ pkg/R/informationCriterion.R	2008-12-02 22:31:58 UTC (rev 69)
@@ -21,7 +21,7 @@
 ## Just returns AIC, AICc, and BIC weights for each of the trees analyzed in a hansenBatch object
   outdata <- vector("list", length(hansenBatch$hansens))
   N = hansenBatch$N
-  for(i in 1:length(outdata)) {
+  for(i in seq(outdata)) {
     temp <- hansenBatch$hansens[[i]]
     outdata[[i]] <- informationCriterion(lnL = temp[, 'loglik'], K = temp[, 'dof'], n = N, names = row.names(temp))
     }

Modified: pkg/R/summarizingAnalyses.R
===================================================================
--- pkg/R/summarizingAnalyses.R	2008-12-02 21:43:17 UTC (rev 68)
+++ pkg/R/summarizingAnalyses.R	2008-12-02 22:31:58 UTC (rev 69)
@@ -32,8 +32,8 @@
     for(i in seq(nnodes)) {
       modelsMatrixSubset <- modelsMatrix[[tree]][hansenBatch$regMatrix$overall[, nodes[i]] == 1, ] # subset models that contain node i
       if(identical(dim(modelsMatrixSubset), NULL)) # is modelsMatrixSubset a 1-d vector? if so then:
-        nodeWeightsSummed[, nodes[i]] <- sum(nodeWeightsSummed[, nodes[i]], modelsMatrixSubset, na.rm = T) # because extracting a single row yields a vector, and dim returns NULL for a vector
-      else nodeWeightsSummed[, nodes[i]] <- sum(nodeWeightsSummed[, nodes[i]], colSums(modelsMatrixSubset, na.rm = T), na.rm = T)
+        nodeWeightsSummed[, nodes[i]] <- nodeWeightsSummed[, nodes[i]] + replace.matrix(modelsMatrixSubset, NA, 0) # because extracting a single row yields a vector, and dim returns NULL for a vector
+      else nodeWeightsSummed[, nodes[i]] <- nodeWeightsSummed[, nodes[i]] + colSums(modelsMatrixSubset, na.rm = T)
 	  }
     sigmaSqVector[tree] <- weighted.mean(hansenBatch$hansens[[tree]][, 'sigma.squared'], bic, na.rm = T)
     if(hansenBatch$brown) bicOU <- bic[1: (length(bic) - 1)]
@@ -49,9 +49,8 @@
                                  )
                                  
   }
+  # in this matrix, the weight for each node is averaged only over trees that possess that node; weights may not sum to 1.0
   nodeWeightsMatrix.unnormalized <- nodeWeightsSummed / matrix(nodeSums, nrow = dim(nodeWeightsSummed)[1], ncol = nnodes, byrow = T)
-                                                          # in this matrix, the weight for each node is averaged only over trees
-                                                          # that possess that node; weights may not sum to 1.0
   nodeWeightsMatrix.allNodes <- nodeWeightsSummed / ntrees 
                                                # in this matrix, the weight for each node is averaged over all trees, setting weight
                                                # equal to zero in any trees that lack that node. Weights sum to 1.0, and they
@@ -95,6 +94,8 @@
   message("ESTIMATED SUPPORT FOR CHANGES OCCURRING AT DESIGNATED NODES")
   message("Averaged over all trees, thus multiplying the clade distribution by the mean support for the model:")
   print(hansenSummary$nodeWeightsMatrix$allNodes)
+  message("\nSupport conditioned on trees that possess the node")
+  print(hansenSummary$nodeWeightsMatrix$unnormalized)
   #message("\nESTIMATED SUPPORT FOR NUMBER OF PARAMETERS IN THE MODEL")
   #message("The properties of this support value have not been studied and are likely to be biased strongly toward the median value of\nK, as K is largest at the median values (they are distributed according to Stirling numbers of the first kind).")
   #print(hansenSummary$kMatrix)