#' @title Linear regression of X,Y,Z-variables with correlated errors
#'
#' @description
#' Implements the maximum likelihood algorithm of Ludwig and
#' Titterington (1994) for linear regression of three dimensional data
#' with correlated uncertainties.
#'
#' @details
#' Ludwig and Titterington (1994)'s 3-dimensional linear regression
#' algorithm for data with correlated uncertainties is an extension of
#' the 2-dimensional algorithm by Titterington and Halliday (1979),
#' which itself is equivalent to the algorithm of York et al. (2004).
#' Given \eqn{n} triplets of (approximately) collinear measurements
#' \eqn{X_i}, \eqn{Y_i} and \eqn{Z_i} (for \eqn{1 \leq i \leq n}),
#' their uncertainties \eqn{s[X_i]}, \eqn{s[Y_i]} and \eqn{s[Z_i]},
#' and their covariances cov[\eqn{X_i,Y_i}], cov[\eqn{X_i,Z_i}] and
#' cov[\eqn{Y_i,Z_i}], the \code{titterington} function fits two
#' slopes and intercepts with their uncertainties. It computes the
#' MSWD as a measure of under/overdispersion. Overdispersed datasets
#' (MSWD>1) can be dealt with in the same three ways that are
#' described in the documentation of the \code{\link{isochron}}
#' function.
#'
#' @param x an \code{[nx9]} matrix with the following columns:
#' \code{X, sX,} \code{Y, sY,} \code{Z, sZ}, \code{rXY,}
#' \code{rXZ,} \code{rYZ}.
#' @return A four-element list of vectors containing:
#'
#' \describe{
#' \item{par}{4-element vector \code{c(a,b,A,B)} where \code{a} is
#' the intercept of the \code{X-Y} regression, \code{b}
#' is the slope of the \code{X-Y} regression, \code{A}
#' is the intercept of the \code{X-Z} regression, and
#' \code{B} is the slope of the \code{X-Z} regression.}
#'
#' \item{cov}{\code{[4x4]}-element covariance matrix of \code{par}}
#'
#' \item{mswd}{the mean square of the residuals (a.k.a `reduced
#' Chi-square') statistic}
#'
#' \item{p.value}{p-value of a Chi-square test for linearity}
#'
#' \item{df}{the number of degrees of freedom for the
#' Chi-square test (2\eqn{n}-4)}
#'
#' \item{tfact}{the \eqn{100(1-\alpha/2)\%} percentile of the
#' t-distribution with \eqn{(n-2k+1)} degrees of freedom}
#' }
#' @examples
#' d <- matrix(c(0.1677,0.0047,1.105,0.014,0.782,0.015,0.24,0.51,0.33,
#' 0.2820,0.0064,1.081,0.013,0.798,0.015,0.26,0.63,0.32,
#' 0.3699,0.0076,1.038,0.011,0.819,0.015,0.27,0.69,0.30,
#' 0.4473,0.0087,1.051,0.011,0.812,0.015,0.27,0.73,0.30,
#' 0.5065,0.0095,1.049,0.010,0.842,0.015,0.27,0.76,0.29,
#' 0.5520,0.0100,1.039,0.010,0.862,0.015,0.27,0.78,0.28),
#' nrow=6,ncol=9)
#' colnames(d) <- c('X','sX','Y','sY','Z','sZ','rXY','rXZ','rYZ')
#' titterington(d)
#' @seealso \code{\link{york}}, \code{\link{isochron}}, \code{\link{ludwig}}
#' @references
#' Ludwig, K.R. and Titterington, D.M., 1994. Calculation
#' of \eqn{^{230}}Th/U isochrons, ages, and errors. Geochimica et
#' Cosmochimica Acta, 58(22), pp.5031-5042.
#'
#' Titterington, D.M. and Halliday, A.N., 1979. On the fitting of
#' parallel isochrons and the method of maximum likelihood. Chemical
#' Geology, 26(3), pp.183-195.
#'
#' York, D., Evensen, N.M., Martinez, M.L. and De Basebe Delgado, J., 2004.
#' Unified equations for the slope, intercept, and standard
#' errors of the best straight line. American Journal of Physics,
#' 72(3), pp.367-375.
#' @export
titterington <- function(x){
ns <- nrow(x)
fitXY <- york(subset(x,select=c(1,2,3,4,7)))
a <- fitXY$a[1]
b <- fitXY$b[1]
fitXZ <- york(subset(x,select=c(1,2,5,6,8)))
A <- fitXZ$a[1]
B <- fitXZ$b[1]
init <- c(a,b,A,B)
dat <- matrix2covlist(x)
fit <- stats::optim(init,fn=L_tit,method="BFGS",gr=gr_tit,dat=dat)
fish <- fisher_tit(fit$par,dat)
AA <- fish[1:ns,1:ns]
BB <- fish[1:ns,(ns+1):(ns+4)]
CC <- fish[(ns+1):(ns+4),1:ns]
DD <- fish[(ns+1):(ns+4),(ns+1):(ns+4)]
out <- list()
out$par <- fit$par
out$cov <- blockinverse(AA,BB,CC,DD)
mswd <- mswd_tit(fit$par,dat)
out <- c(out,mswd)
parnames <- c('a','b','A','B')
names(out$par) <- parnames
rownames(out$cov) <- parnames
colnames(out$cov) <- parnames
out$type <- 'titterington'
out
}
get_titterington_xyz <- function(XYZ,a,b,A,B,w=0,wtype=NA){
vX <- XYZ[,'sX']^2
vY <- XYZ[,'sY']^2
vZ <- XYZ[,'sZ']^2
if (wtype=='a') vX <- vX + w^2
if (wtype=='A') vZ <- vZ + w^2
sXY <- XYZ[,'rXY']*XYZ[,'sX']*XYZ[,'sY']
sXZ <- XYZ[,'rXZ']*XYZ[,'sX']*XYZ[,'sZ']
sYZ <- XYZ[,'rYZ']*XYZ[,'sY']*XYZ[,'sZ']
omega <- invertcovmat(vx=vX,vy=vY,vz=vZ,sxy=sXY,sxz=sXZ,syz=sYZ)
r <- XYZ[,'Y']-a-b*XYZ[,'X']
R <- XYZ[,'Z']-A-B*XYZ[,'X']
alpha <- omega[,'xx'] + omega[,'yy']*b^2 + omega[,'zz']*B^2 +
2*omega[,'xy']*b + 2*omega[,'xz']*B + 2*omega[,'yz']*b*B
beta <- r*(omega[,'xy'] + b*omega[,'yy'] + B*omega[,'yz']) +
R*(omega[,'xz'] + b*omega[,'yz'] + B*omega[,'zz'])
if (wtype%in%c('b','B')){
slopetitroot <- function(p,XYZi,a,b,A,B,w){
E <- dEdx <- matrix(0,3,3)
E[1,1] <- XYZi['sX']^2
if (wtype=='b'){
E[2,2] <- XYZi['sY']^2 + (w*p)^2
dEdx[2,2] <- 2*p*w^2
} else {
E[2,2] <- XYZi['sY']^2
}
if (wtype=='B') {
E[3,3] <- XYZi['sZ']^2 + (w*p)^2
dEdx[3,3] <- 2*p*w^2
} else {
E[3,3] <- XYZi['sZ']^2
}
dDdx <- c(-1,-b,-B)
E[1,2] <- E[2,1] <- XYZi['rXY']*XYZi['sX']*XYZi['sY']
E[1,3] <- E[3,1] <- XYZi['rXZ']*XYZi['sX']*XYZi['sZ']
E[2,3] <- E[3,2] <- XYZi['rYZ']*XYZi['sY']*XYZi['sZ']
D <- c(XYZi['X']-p,XYZi['Y']-a-b*p,XYZi['Z']-A-B*p)
O <- solve(E)
sum(diag(O%*%dEdx)) + 2*D%*%O%*%dDdx - D%*%O%*%dEdx%*%O%*%D
}
slopetithelper <- function(i,XYZ,a,b,A,B,w,lower,upper){
stats::uniroot(f=slopetitroot,lower=lower[i],upper=upper[i],
XYZi=XYZ[i,],a=a,b=b,A=A,B=B,w=w,extendInt='yes',
tol=.Machine$double.eps^0.5)
}
ns <- nrow(XYZ)
init <- XYZ[,'X'] + beta/alpha
dx <- diff(range(XYZ[,'X']))
fit <- sapply(1:ns,slopetithelper,lower=init-dx,
upper=init+dx,XYZ=XYZ,a=a,b=b,A=A,B=B,w=w)
x <- unlist(fit[1,])
} else {
x <- XYZ[,'X'] + beta/alpha
}
y <- a+b*x
z <- A+B*x
cbind('x'=x,'y'=y,'z'=z)
}
mswd_tit <- function(abAB,dat){
a <- abAB[1]
b <- abAB[2]
A <- abAB[3]
B <- abAB[4]
XYZ <- dat$XYZ
omega <- dat$omega
ns <- length(XYZ)
S <- 0
for (i in 1:ns){
X <- XYZ[[i]][1]
Y <- XYZ[[i]][2]
Z <- XYZ[[i]][3]
abg <- alpha_beta_gamma(a,b,A,B,XYZ[[i]],omega[[i]])
alpha <- abg[1]
beta <- abg[2]
gamma <- abg[3]
x <- X + beta/alpha
S <- S + alpha*(X-x)^2 + 2*beta*(X-x) + gamma
}
out <- list()
out$df <- 2*ns-4
out <- append(out,getMSWD(S,out$df))
out
}
matrix2covlist <- function(x){
ns <- nrow(x)
out <- list()
out$XYZ <- list()
out$omega <- list()
covmat <- matrix(0,3,3)
for (i in 1:ns){
covmat <- cor2cov3(x[i,2],x[i,4],x[i,6],x[i,7],x[i,8],x[i,9])
out$XYZ[[i]] <- x[i,c(1,3,5)]
out$omega[[i]] <- solve(covmat)
}
out
}
alpha_beta_gamma <- function(a,b,A,B,XYZ,omega){
r <- XYZ[2]-a-b*XYZ[1]
R <- XYZ[3]-A-B*XYZ[1]
alpha <- omega[1,1] + omega[2,2]*b^2 + omega[3,3]*B^2 +
2*omega[1,2]*b + 2*omega[1,3]*B + 2*omega[2,3]*b*B
beta <- r*(omega[1,2] + b*omega[2,2] + B*omega[2,3]) +
R*(omega[1,3] + b*omega[2,3] + B*omega[3,3])
gamma <- omega[2,2]*r^2 + omega[3,3]*R^2 + 2*omega[2,3]*r*R
c(alpha,beta,gamma)
}
fisher_tit <- function(abAB,dat){
a <- abAB[1] # y intercept
b <- abAB[2] # y slope
A <- abAB[3] # z intercept
B <- abAB[4] # z slope
ns <- length(dat$XYZ)
out <- matrix(0,ns+4,ns+4)
d2L.da2 <- 0
d2L.dadb <- 0
d2L.db2 <- 0
d2L.dA2 <- 0
d2L.dAdB <- 0
d2L.dB2 <- 0
d2L.dadA <- 0
d2L.dadB <- 0
d2L.dbdA <- 0
d2L.dbdB <- 0
for (i in 1:ns){
XYZ <- dat$XYZ[[i]]
X <- XYZ[1]
Y <- XYZ[2]
Z <- XYZ[3]
omega <- dat$omega[[i]]
abg <- alpha_beta_gamma(a,b,A,B,XYZ,omega)
alpha <- abg[1]
beta <- abg[2]
gamma <- abg[3]
x <- X + beta/alpha
d2L.da2 <- d2L.da2 - omega[2,2]
d2L.dadb <- d2L.dadb - x*omega[2,2]
d2L.db2 <- d2L.db2 - omega[2,2]*x^2
d2L.dA2 <- d2L.dA2 - omega[3,3]
d2L.dAdB <- d2L.dAdB - x*omega[3,3]
d2L.dB2 <- d2L.dB2 - omega[3,3]*x^2
d2L.dadA <- d2L.dadA - omega[2,3]
d2L.dadB <- d2L.dadB - x*omega[2,3]
d2L.dbdA <- d2L.dbdA - x*omega[2,3]
d2L.dbdB <- d2L.dbdB - omega[2,3]*x^2
d2L.dadx <- -(omega[1,2] + b*omega[2,2] + B*omega[2,3])
d2L.dbdx <- -x*(omega[1,2] + b*omega[2,2] + B*omega[2,3])
d2L.dAdx <- -(omega[1,3] + b*omega[2,3] + B*omega[3,3])
d2L.dBdx <- -x*(omega[1,3] + b*omega[2,3] + B*omega[3,3])
d2L.dx2 <- -(omega[1,1] + 2*b*omega[1,2] + 2*B*omega[1,3] +
omega[2,2]*b^2 + omega[2,3]*B^2 + 2*omega[3,3]*b*B)
out[i,i] <- -d2L.dx2
out[i,ns+1] <- -d2L.dadx
out[i,ns+2] <- -d2L.dbdx
out[i,ns+3] <- -d2L.dAdx
out[i,ns+4] <- -d2L.dBdx
out[ns+1,i] <- out[i,ns+1]
out[ns+2,i] <- out[i,ns+2]
out[ns+3,i] <- out[i,ns+3]
out[ns+4,i] <- out[i,ns+4]
}
out[ns+1,ns+1] <- -d2L.da2
out[ns+2,ns+2] <- -d2L.db2
out[ns+3,ns+3] <- -d2L.dA2
out[ns+4,ns+4] <- -d2L.dB2
out[ns+1,ns+2] <- -d2L.dadb
out[ns+1,ns+3] <- -d2L.dadA
out[ns+1,ns+4] <- -d2L.dadB
out[ns+2,ns+3] <- -d2L.dbdA
out[ns+2,ns+4] <- -d2L.dbdB
out[ns+3,ns+4] <- -d2L.dAdB
out[ns+2,ns+1] <- out[ns+1,ns+2]
out[ns+3,ns+1] <- out[ns+1,ns+3]
out[ns+3,ns+2] <- out[ns+2,ns+3]
out[ns+4,ns+1] <- out[ns+1,ns+4]
out[ns+4,ns+2] <- out[ns+2,ns+4]
out[ns+4,ns+3] <- out[ns+3,ns+4]
out
}
gr_tit <- function(abAB,dat){
ns <- length(dat)
a <- abAB[1] # y intercept
b <- abAB[2] # y slope
A <- abAB[3] # z intercept
B <- abAB[4] # z slope
dS.da <- 0
dS.db <- 0
dS.dA <- 0
dS.dB <- 0
dL.dS <- -1/2
dr.da <- -1
dR.dA <- -1
dR.da <- 0
dr.dA <- 0
dR.db <- 0
dr.dB <- 0
dalpha.da <- 0
dalpha.dA <- 0
for (i in 1:ns){
XYZ <- dat$XYZ[[i]]
X <- XYZ[1]
Y <- XYZ[2]
Z <- XYZ[3]
omega <- dat$omega[[i]]
r <- Y-a-b*X
R <- Z-A-B*X
abg <- alpha_beta_gamma(a,b,A,B,XYZ,omega)
alpha <- abg[1]
beta <- abg[2]
gamma <- abg[3]
dr.db <- -X
dR.dB <- -X
dbeta.da <- dr.da*(omega[1,2] + b*omega[2,2] + B*omega[2,3]) +
dR.da*(omega[1,3] + b*omega[2,3] + B*omega[3,3])
dgamma.da <- 2*(omega[2,2]*r*dr.da + omega[3,3]*R*dR.da +
omega[2,3]*R*dr.da + omega[2,3]*r*dR.da)
dalpha.db <- 2*(omega[2,2]*b + omega[1,2] + omega[2,3]*B)
dbeta.db <- dr.db*(omega[1,2] + b*omega[2,2] + B*omega[2,3]) +
dR.db*(omega[1,3] + b*omega[2,3] + B*omega[3,3]) +
r*omega[2,2] + R*omega[2,2]
dgamma.db <- 2*(omega[2,2]*r*dr.db + omega[3,3]*R*dR.db +
omega[2,3]*R*dr.db + omega[2,3]*r*dR.db)
dbeta.dA <- dr.dA*(omega[1,2] + b*omega[2,2] + B*omega[2,3]) +
dR.dA*(omega[1,3] + b*omega[2,3] + B*omega[3,3])
dgamma.dA <- 2*(omega[2,2]*r*dr.dA + omega[3,3]*R*dR.dA +
omega[2,3]*R*dr.dA + omega[2,3]*r*dR.dA)
dalpha.dB <- 2*(omega[3,3]*B + omega[1,3] + omega[2,3]*b)
dbeta.dB <- dr.dB*(omega[1,2] + b*omega[2,2] + B*omega[2,3]) +
dR.dB*(omega[1,3] + b*omega[2,3] + B*omega[3,3]) +
r*omega[2,3] + R*omega[3,3]
dgamma.dB <- 2*(omega[2,2]*r*dr.dB + omega[3,3]*R*dR.dB +
omega[2,3]*R*dr.dB + omega[2,3]*r*dR.dB)
dS.da <- dgamma.da - 2*(beta/alpha)*dbeta.da + dalpha.da*(beta/alpha)^2
dS.db <- dgamma.db - 2*(beta/alpha)*dbeta.db + dalpha.db*(beta/alpha)^2
dS.dA <- dgamma.dA - 2*(beta/alpha)*dbeta.dA + dalpha.dA*(beta/alpha)^2
dS.dB <- dgamma.dB - 2*(beta/alpha)*dbeta.dB + dalpha.dB*(beta/alpha)^2
}
c(dS.da,dS.db,dS.dA,dS.dB)*dL.dS
}
L_tit <- function(abAB,dat){
ns <- length(dat)
a <- abAB[1] # y intercept
b <- abAB[2] # y slope
A <- abAB[3] # z intercept
B <- abAB[4] # z slope
S <- 0 # initialise sum of squares
for (i in 1:ns){
abg <- alpha_beta_gamma(a,b,A,B,dat$XYZ[[i]],dat$omega[[i]])
alpha <- abg[1]
beta <- abg[2]
gamma <- abg[3]
S <- S + gamma - (beta^2)/alpha
}
S/2
}
#' Convert IsoplotR data to Titterington ratios
#' @noRd
data2tit <- function(x,...){
UseMethod("data2tit",x)
}
#' @noRd
data2tit.default <- function(x,...){
stop('default function undefined')
}
#' @param osmond if \code{TRUE}, generates 2/8 - 4/8 - 0/8 ratios
#' if \code{FALSE}, generates 8/2 - 4/2 - 0/2 ratios
#' @param generic \if \code{TRUE}, uses X, Y, Z labels,
#' otherwise uses isotopic ratio labels
#' @noRd
data2tit.ThU <- function(x,osmond=TRUE,generic=TRUE,...){
ns <- length(x)
out <- matrix(0,ns,9)
if ((x$format == 1 & !osmond) | (x$format == 2 & osmond)){
out <- x$x
} else {
for (i in 1:ns){
out[i,] <- ThU_convert(x$x[i,])
}
}
if (generic){
colnames(out) <- c('X','sX','Y','sY','Z','sZ','rXY','rXZ','rYZ')
} else if (osmond){
colnames(out) <- c('Th232U238','sTh232U238','U234U238','sU234U238',
'Th230U238','sTh230U238','rXY','rXZ','rYZ')
} else {
colnames(out) <- c('U238Th232','sU238Th232','U234Th232','su234Th232',
'Th230Th232','sTh230Th232','rXY','rXZ','rYZ')
}
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.