Dear R-helpers,
 
I´m in the context of writing a general function for error propagation
in R.
There are somehow a few questions I would like to ask (discuss), as my
statistical knowledge is somewhat restricted.
Below is the function I wrote, the questions are marked.
Many thanks in advance.
 
propagate <- function(expr, varList, type = c("stat", "raw"), cov =
TRUE)
{
      require(gtools, quietly = TRUE)
      if (!is.expression(expr)) stop("not a valid expression!")
      if (!is.list(varList)) stop("values are not of type 'list'!")
      type <- match.arg(type)
 
      m <- match(all.vars(expr), names(varList))
      if (any(is.na(m))) stop(paste("some variables do not match"))
      lenCheck <- sapply(varList, function(x) length(x))
      if (sum(lenCheck == 2) == length(varList) && type == "raw")
print("All list items of size 2. Are you sure this is raw data?")
      
      varList <- varList[m]
      BT <- NULL
      COVALL <- 0
 
      ### Q1: These are the permutations for the covariance matrix. Is
this right? No repeats allowed because cov(X1, X1) equals variance?
      perm <- permutations(length(varList), 2, repeats.allowed = FALSE)
 
      ### Q2: A Bartletts test to check for hemogeneity of variance,
which is prerequisite for gaussian error propagation. Does that make
sense? 
      if (type == "raw") {
            BT <- bartlett.test(varList)
            if (BT$p.value < 0.05) print("Bartlett's test indicates
heteroscedasticity! Continuing anyway, but keep in mind...")
            means <- sapply(varList, function(x)  mean(x, na.rm = TRUE))
            sds <- sapply(varList, function(x)  sd(x, na.rm = TRUE))
            VARS <- as.data.frame(rbind(means, sds))
      }
      else {
            means <- sapply(varList, function(x) x[1])
            sds <- sapply(varList, function(x)  x[2])
            VARS <- as.data.frame(rbind(means, sds))
      }
 
     ### partial derivatives for each variable in function
      el <- lapply(all.vars(expr), function(x) deriv(expr, x))
 
      part <- vector()
      term <- vector()
 
     ### evaluate partial derivatives with the values from the
variables.
      for (i in 1:length(el)) {
            part[i] <- attr(eval(el[[i]], envir = VARS[1, ]),
"gradient")
            term[i] <- (part[i] * VARS[2, i])^2
      }
                                           
      ### Q3: This is the part I´m not sure about. Often the covariance
matrix is neglected in the literature, but I read this is not good
practice.
     ### The covariances are evaluated for each pair of the permutation,
but without pairs with the same index.
     ### Is this right or am I m issing something here???
      if (cov) {
            for (i in 1:nrow(perm)) {
                  COV01 <- part[perm[i, 1]]
                  COV02 <- part[perm[i, 2]]
                  COV03 <- as.vector(varList[[perm[i, 1]]])
                  COV04 <- as.vector(varList[[perm[i, 2]]])
                  COVALL[i] <- COV01 * COV02 * cov(COV03, COV04, use =
"complete.obs")
            }
      }
      
      covMat <- cov(as.data.frame(varList))
      ### calculate error from square root af all squared terms
      ERROR <- sqrt(sum(term) + sum(COVALL))
      return(list(error = ERROR, partials = el, htest = BT, covMat =
covMat))
}
 
############
Example:
e <- expression((E1^cp1)/(E2^cp2))
data <- list(E1 = c(1.8, 1.7, 1.8), E2 = c(1.65, 1.72, 1.71), cp1 =
c(17.6, 17.8, 17.7), cp2 = c(22.12, 21.5, 22.1))
temp <- propagate(e, data, type = "raw", cov = T)
print(temp)
 
Many thanks in advance!!
 
 
 
_____________________________________
Dr. rer. nat. Andrej-Nikolai Spiess (Dipl. Biol.)
Department of Andrology
University Hospital Hamburg-Eppendorf
Martinistr. 52
20246 Hamburg
phone: +49-40-428031585
email: a.spiess@uke.uni-hamburg.de
 
 


-- 
Pflichtangaben gemäß Gesetz über elektronische Handelsregister und Genossenschaftsregister sowie das Unternehmensregister (EHUG):

Universitätsklinikum Hamburg-Eppendorf
Körperschaft des öffentlichen Rechts
Gerichtsstand: Hamburg

Vorstandsmitglieder:
Prof. Dr. Jörg F. Debatin (Vorsitzender)
Dr. Alexander Kirstein
Ricarda Klein
Prof. Dr. Dr. Uwe Koch-Gromus

	[[alternative HTML version deleted]]

