r-gregmisc-users

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

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

Revision: 1525
          http://r-gregmisc.svn.sourceforge.net/r-gregmisc/?rev=1525&view=rev
Author:   warnes
Date:     2012-04-19 22:05:48 +0000 (Thu, 19 Apr 2012)
Log Message:
-----------
More fixes for support of S4 'mer' class from lme4 package.

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

Modified: trunk/gmodels/R/ci.R
===================================================================
--- trunk/gmodels/R/ci.R	2012-04-19 21:13:28 UTC (rev 1524)
+++ trunk/gmodels/R/ci.R	2012-04-19 22:05:48 UTC (rev 1525)
@@ -85,7 +85,10 @@
                                 show.beta0 = FALSE,
                                 n.sim = n.sim)
 
-    retval <- retval[, c("Estimate", "Lower.CI", "Upper.CI", "Std. Error", "p value")]
+    retval <- retval[,
+                     c("Estimate", "Lower.CI", "Upper.CI", "Std. Error", "p value"),
+                     drop=FALSE
+                     ]
     colnames(retval)[c(2:3, 5)] <- c("CI lower", "CI upper", "p-value")
     rownames(retval) <- names(x.effects)
     

Modified: trunk/gmodels/R/est.mer.R
===================================================================
--- trunk/gmodels/R/est.mer.R	2012-04-19 21:13:28 UTC (rev 1524)
+++ trunk/gmodels/R/est.mer.R	2012-04-19 22:05:48 UTC (rev 1525)
@@ -4,6 +4,7 @@
 # Laboratory of Genomic Diversity at NCI Frederick
 # SAIC Frederick, Inc
 # Created April 25, 2006
+# Updated 2012-04-19 for S4 version of lmer object
 
 est.mer <- function(obj, cm, beta0, conf.int, show.beta0, n.sim)
 {
@@ -17,32 +18,34 @@
     n <- dim(cm)[2]
 
   ## drop extra information on end
-  samp.cm <- as.matrix(samp@fixef)[, 1:n] %*% t(cm)
+  samp.cm <- t(samp@fixef)[1:n,, drop=FALSE] %*% cm
 
   # calculate requested statistics
   est <- apply(samp.cm, 2, mean)
-  stderr <- sd(samp.cm)
+  stderr <- apply(samp.cm, 2, sd)
   
   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)
+    {
+      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,]
-  }
+      lower.ci <- samp.ci[1,]
+      upper.ci <- samp.ci[2,]
+    }
 
   # return results
   if(!show.beta0)
@@ -61,8 +64,11 @@
   return(samp.stats)
 }
 
-percentile <- function(x, distn)
+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.

[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.