[R] Error propagation

From: Andrej-Nikolai Spiess <a.spiess_at_uke.uni-hamburg.de>
Date: Tue, 25 Mar 2008 13:30:20 +0100


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_at_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]]

______________________________________________ R-help_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.

Received on Tue 25 Mar 2008 - 13:15:30 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 25 Mar 2008 - 14:30:23 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive