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: [EMAIL PROTECTED]
--
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@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.
