# [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)
sds <- sapply(varList, function(x)  x)
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, ]),
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)

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