[R-gregmisc-users] SF.net SVN: r-gregmisc:[1334] trunk/gmodels/R/est.mer.R

[<wa...@us...>]

Revision: 1334
          http://r-gregmisc.svn.sourceforge.net/r-gregmisc/?rev=1334&view=rev
Author:   warnes
Date:     2009-05-09 05:01:47 +0000 (Sat, 09 May 2009)

Log Message:
-----------
Add support for lme4's 'mer' objects

Added Paths:
-----------
    trunk/gmodels/R/est.mer.R

Added: trunk/gmodels/R/est.mer.R
===================================================================
--- trunk/gmodels/R/est.mer.R	                        (rev 0)
+++ trunk/gmodels/R/est.mer.R	2009-05-09 05:01:47 UTC (rev 1334)
@@ -0,0 +1,70 @@
+# est.mer.R
+# generate estimable output for mer objects using mcmcsamp()
+# Randall Johnson
+# Laboratory of Genomic Diversity at NCI Frederick
+# SAIC Frederick, Inc
+# Created April 25, 2006
+
+est.mer <- function(obj, cm, beta0, conf.int, show.beta0, n.sim)
+{
+##   if(!require(coda, quietly=TRUE))
+##     stop("coda package required when sim.mer == TRUE")
+
+  samp <- mcmcsamp(obj, n.sim)
+##  samp.summ <- summary(samp)
+
+  if(is.null(dim(cm)))
+    n <- length(cm)
+  else
+    n <- dim(cm)[2]
+                               # drop extra information on end
+  samp.cm <- as.matrix(samp)[, 1:n] %*% t(cm)
+
+  # calculate requested statistics
+  est <- apply(samp.cm, 2, mean)
+  stderr <- sd(samp.cm)
+  
+  pval <- sapply(1:length(beta0),
+                 function(i){percentile(beta0[i], samp.cm[,i])})
+  pval <- ifelse(pval <= .5, 2*pval, 2*(1-pval))
+
+  if(is.null(conf.int))
+  {
+    lower.ci <- NULL
+    upper.ci <- NULL
+  }else{
+    alpha <- 1-conf.int
+    samp.ci <- sapply(1:length(beta0),
+                      function(i)
+                        {
+                          quantile(samp.cm[,i], probs=c(alpha/2, 1-alpha/2))
+                        }
+                      )
+
+    lower.ci <- samp.ci[1,]
+    upper.ci <- samp.ci[2,]
+  }
+
+  # return results
+  if(!show.beta0)
+    beta0 <- NULL
+  
+  samp.stats <- cbind('beta0' = beta0,
+                      'Estimate' = est,
+                      'Std. Error' = stderr,
+                      'p value' = pval,
+                      'Lower.CI' = lower.ci,
+                      'Upper.CI' = upper.ci)
+
+  row.names(samp.stats) <- paste('(', apply(cm, 1, paste, collapse=" "),
+                                 ')', sep='')
+    
+  return(samp.stats)
+}
+
+percentile <- function(x, distn)
+{
+  n <- length(distn)
+    
+  return(findInterval(x, distn[order(distn)]) / n)
+}


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.