[R-gregmisc-users] SF.net SVN: r-gregmisc:[1575] trunk/gmodels/R

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

Revision: 1575
          http://r-gregmisc.svn.sourceforge.net/r-gregmisc/?rev=1575&view=rev
Author:   warnes
Date:     2012-06-28 00:41:26 +0000 (Thu, 28 Jun 2012)
Log Message:
-----------
Update est.mer() to support new S4 "mer" class.

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

Removed Paths:
-------------
    trunk/gmodels/R/est.lmer.R

Deleted: trunk/gmodels/R/est.lmer.R
===================================================================
--- trunk/gmodels/R/est.lmer.R	2012-06-28 00:40:29 UTC (rev 1574)
+++ trunk/gmodels/R/est.lmer.R	2012-06-28 00:41:26 UTC (rev 1575)
@@ -1,70 +0,0 @@
-# est.lmer.R
-# generate estimable output for lmer objects using mcmcsamp()
-# Randall Johnson
-# Laboratory of Genomic Diversity at NCI Frederick
-# SAIC Frederick, Inc
-# Created April 25, 2006
-
-est.lmer <- function(obj, cm, beta0, conf.int, show.beta0, n.sim)
-{
-##  if(!require(coda, quietly=TRUE))
-##    stop("coda package required when sim.lmer == 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 <- drop(cm %*% samp.summ$statistics[1:n,1])
-  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)
-}

Modified: trunk/gmodels/R/est.mer.R
===================================================================
--- trunk/gmodels/R/est.mer.R	2012-06-28 00:40:29 UTC (rev 1574)
+++ trunk/gmodels/R/est.mer.R	2012-06-28 00:41:26 UTC (rev 1575)
@@ -9,17 +9,11 @@
 est.mer <- function(obj, cm, beta0, conf.int, show.beta0, n.sim)
 {
 
-  samp <- mcmcsamp(obj, n.sim)
+  samp <- lme4:::mcmcsamp(obj, n.sim)
   ##  samp.summ <- summary(samp)
 
-  if(is.null(dim(cm)))
-    n <- length(cm)
-  else
-    n <- dim(cm)[2]
+  samp.cm <- t(cm %*% samp@fixef)
 
-  ## drop extra information on end
-  samp.cm <- t(samp@fixef)[1:n,, drop=FALSE] %*% cm
-
   # calculate requested statistics
   est <- apply(samp.cm, 2, mean)
   stderr <- apply(samp.cm, 2, sd)
@@ -64,12 +58,3 @@
   return(samp.stats)
 }
 
-percentile <- function(x, distn, include.observed=FALSE)
-{
-  if(include.observed)
-    distn <- c(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.