Nothing
#' Extrapolation to 'time zero'
#'
#' This function extrapolates time resolved mass spectrometer data to
#' t=0. When fed with multicollector data, it forms the ratios of the
#' raw signals, forms their logs and performs linear regression to t=0
#' When fed with single collector data, the function first takes their
#' logs and extrapolates them to t=0 before taking ratios, unless
#' \code{denmass}=NULL, in which case the logs of the raw signals are
#' extrapolated.
#'
#' @param x an object of class \code{timeresolved} or \code{PHdata}
#' @param ... further arguments (see below)
#' @return an object of class \code{logratios}
#' @examples
#' samplefile <- system.file("Samples.csv",package="ArArRedux")
#' masses <- c("Ar37","Ar38","Ar39","Ar40","Ar36")
#' m <- loaddata(samplefile,masses) # samples and J-standards
#' blanklabel <- "EXB#"
#' l <- fitlogratios(blankcorr(m,blanklabel),"Ar40")
#' plotcorr(l)
#' @export
fitlogratios <- function(x,...){ UseMethod("fitlogratios",x) }
#' @rdname fitlogratios
#' @export
fitlogratios.default <- function(x,...){ stop() }
#' @param denmass a string denoting the denominator isotope
#' @rdname fitlogratios
#' @export
fitlogratios.timeresolved <- function(x,denmass,...){
r <- takeratios(x,denmass)
l <- takelogs(r)
f <- timezero(l)
return(cast(f,"logratios"))
}
#' @rdname fitlogratios
#' @export
fitlogratios.PHdata <- function(x,denmass=NULL,...){
for (i in 1:nmasses(x)){
z <- newfit(x)
l <- takelogs(x$signals[[i]])
z <- setmasses(z,z$masses[i],timezero(l))
}
if (is.null(denmass)) {
f <- z
f$denmass <- NA
} else {
f <- takeratios(z,denmass)
}
return(cast(f,"logratios"))
}
Jtakeratios <- function(nruns,inum,iden){
nlr <- length(inum)
nmasses <- nlr+1
J <- matrix(0,nrow=nlr,ncol=nmasses) # elementary matrix
for (i in 1:length(inum)){
J[i,inum[i]] <- 1
J[,iden] <- -1
}
out <- matrix(0,nrow=nlr*nruns,ncol=nmasses*nruns)
for (i in 1:nruns){
irow <- ((i-1)*nlr+1):(i*nlr)
icol <- ((i-1)*nmasses+1):(i*nmasses)
out[irow,icol] <- J
}
return(out)
}
takeratios <- function(x,...){ UseMethod("takeratios",x) }
takeratios.default <- function(x,...){stop()}
takeratios.timeresolved <- function(x,denmass,...){
den <- getmasses(x,denmass)
num <- getmasses(x,denmass,invert=TRUE)
out <- num
for (mass in num$masses){ # loop through the isotopes
ratio <- getmasses(num,mass)$d/den$d
out <- setmasses(out,mass,ratio)
}
out$denmass <- denmass
class(out) <- append(class(out),"ratio")
return(out)
}
takeratios.fit <- function(x,denmass,...){
out <- x
iden <- which(x$mass == denmass)
inum <- which(x$mass != denmass)
J <- Jtakeratios(nruns(x),inum,iden)
out$masses <- x$masses[inum]
out$intercepts <- J %*% x$intercepts
out$covmat <- J %*% x$covmat %*% t(J)
out$denmass <- denmass
return(out)
}
takelogs <- function(x){
out <- x
out <- replacenegatives(x)
out$d <- log(out$d)
class(out) <- append(class(out),"logged")
return(out)
}
newfit <- function(x,...){ UseMethod("newfit",x) }
newfit.default <- function(x,...){ stop() }
newfit.timeresolved <- function(x,nmasses=NULL,nruns=NULL,...){
out <- x
class(out) <- "fit"
if (is.null(nmasses)) nmasses <- nmasses(x)
if (is.null(nruns)) nruns <- nruns(x)
out$d <- NULL
out$thetime <- NULL
out$intercepts <- rep(0,nmasses*nruns)
out$covmat <- matrix(0,nrow=nmasses*nruns,ncol=nmasses*nruns)
return(out)
}
newfit.PHdata <- function(x,...){
x1 <- x$signals[[1]]
out <- newfit(x1,nmasses(x),nruns(x1))
out$masses <- x$masses
out$signals <- NULL
return(out)
}
timezero <- function(x){
nmasses <- nmasses(x)
out <- newfit(x,nmasses)
bi <- rle(as.vector(x$blankindices))$values # blank indices
irunsout <- 0
for (i in bi){ # loop through the groups
irunsx <- which(i==x$blankindices)
irunsout <- utils::tail(irunsout,n=1) + (1:length(irunsx))
g <- subset(x,irunsx) # extract group
out <- setfit(out,fit(g),nmasses,irunsout)
}
return(out)
}
# x and f are both objects of class "fit"
# nmasses = number of masses or isotope ratios per sample
# iruns = indicates which samples need replacing
setfit <- function(x,f,nmasses,iruns){
out <- x
ii <- getindices(nmasses=nmasses,iruns=iruns)
out$intercepts[ii] <- f$intercepts
out$covmat[ii,ii] <- f$covmat
return(out)
}
fit <- function(x){
out <- newfit(x)
if (nruns(x)>1) { # average time for all samples in group
thetime <- apply(x$thetime,1,"mean")
} else {
thetime <- x$thetime
}
f <- stats::lm(x$d ~ thetime)
if (nruns(x) == 1 & nmasses(x) == 1) { # only one intercept
out$intercepts <- stats::coef(f)["(Intercept)"]
covcolname <- "(Intercept)"
} else { # an entire row of intercepts
out$intercepts <- stats::coef(f)["(Intercept)",]
covcolname <- ":(Intercept)"
}
myvcov <- stats::vcov(f)
j <- which(colnames(myvcov)==covcolname)
out$covmat <- myvcov[j,j]
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.