[Mattice-commits] r62 - pkg/R

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

Author: andrew_hipp
Date: 2008-12-02 08:57:12 +0100 (Tue, 02 Dec 2008)
New Revision: 62

Modified:
   pkg/R/summarizingAnalyses.R
Log:
summary now returns model-averaged alpha and sigma, over trees and models

Modified: pkg/R/summarizingAnalyses.R
===================================================================
--- pkg/R/summarizingAnalyses.R	2008-12-01 21:15:15 UTC (rev 61)
+++ pkg/R/summarizingAnalyses.R	2008-12-02 07:57:12 UTC (rev 62)
@@ -17,6 +17,8 @@
   icObject <- informationCriterion.hansenBatch(hansenBatch)
   nmodels <- dim(hansenBatch$hansens[[1]])[1]
   ntrees <- length(hansenBatch$hansens)
+  sigmaSqVector <- numeric(ntrees)
+  alphaVector <- numeric(ntrees)
   modelsMatrix <- vector('list', ntrees)
   matrixRows <- c('AIC.weight', 'AICc.weight', 'BIC.weight')
   outMatrix <- matrix(0, nrow = length(icObject[[1]]$AICwi), ncol = length(matrixRows), dimnames = list(icObject[[1]]$names, matrixRows))
@@ -29,7 +31,10 @@
     temp <- modelsMatrix[[tree]]; temp[!is.na(temp)] <- 1
     nmodelsMatrix <- nmodelsMatrix + replace.matrix(temp, NA, 0) # nmodelsMatrix gets a 1 for each model present, a 0 for each model not present
     outMatrix <- outMatrix + replace.matrix(modelsMatrix[[tree]], NA, 0) # replace NAs with 0 so that sum works correctly at the end
-  
+    sigmaSqVector[tree] <- weighted.mean(hansenBatch$hansens[[tree]][, 'sigma.squared'], bicList, na.rm = T)
+    alphaVector[tree] <- ifelse(hansenBatch$brown, 
+                                weighted.mean(hansenBatch$hansens[[tree]][1:(nmodels - 1), 'theta / alpha'], bicList[1: (nmodels - 1)], na.rm = T),
+                                weighted.mean(hansenBatch$hansens[[tree]][ , 'theta / alpha'], bicList, na.rm = T) )
     ## the lines below made the weights on branches ignore the fact that the Brownian motion model was part of the
     ##   model set; however, I've removed them b/c support for the Brownian motion model does (and should) contribute 
     ##   to reduced probability of change at any of the nodes. You can uncomment them if you feel differently.
@@ -77,8 +82,10 @@
     if(identical(dim(modelsMatrixSubset), NULL)) kMatrix[, i] <- modelsMatrixSubset # is modelsMatrixSubset a 1-d vector?
     else kMatrix[, i] <- apply(modelsMatrixSubset, 2, sum) 
   }
+  modelAvgAlpha <- mean(alphaVector, na.rm = T)
+  modelAvgSigmaSq <- mean(sigmaSqVector, na.rm = T)
   warning('the K-matrix as currently implemented is calculated over the all-nodes (normalized) weights.')
-  outdata <- list(modelsMatrix = modelsMatrix, weightsMatrix = list(unnormalized = weightsMatrix.unnormalized, allNodes = weightsMatrix.allNodes), nodeWeightsMatrix = list(unnormalized = nodeWeightsMatrix.unnormalized, allNodes = nodeWeightsMatrix.allNodes), kMatrix = kMatrix)
+  outdata <- list(modelsMatrix = modelsMatrix, weightsMatrix = list(unnormalized = weightsMatrix.unnormalized, allNodes = weightsMatrix.allNodes), nodeWeightsMatrix = list(unnormalized = nodeWeightsMatrix.unnormalized, allNodes = nodeWeightsMatrix.allNodes), kMatrix = kMatrix, modelAvgAlpha = modelAvgAlpha, modelAvgSigmaSq = modelAvgSigmaSq)
   class(outdata) <- 'hansenSummary'
   return(outdata)
 }
@@ -99,4 +106,7 @@
   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)
+  message("\nMODEL AVERAGED PARAMETERS")
+  message(paste("alpha =", hansenSummary$modelAvgAlpha))
+  message(paste("sigma^2 =", hansenSummary$modelAvgSigmaSq))
 }
\ No newline at end of file