Nothing
R/agd.r
In AGD: Analysis of Growth Data
Defines functions
adjust.tail
check.names
y2z
z2y
Documented in
y2z
z2y
# agd.r
#
# Tools for course Analysis of Growth Data
### SvB 20dec2014
#'Convert standard deviation scores (SDS) to measurements
#'
#'Converts standard deviation score (SDS) into measurements using
#'an age- and sex-conditional external reference.
#'
#'Functions \code{z2y()} and \code{y2z()} are the inverse of each other.
#'
#'The argument \code{dist} determines the statistical distribution. The
#'possibilities are as follows: \describe{ \item{list("\"NO\"")}{\code{ref}
#'should contain columns \code{mean} and \code{sd}, containing the mean and the
#'standard deviation in the external reference population.}
#'\item{list("\"LMS\"")}{\code{ref} should contain columns \code{L}, \code{S}
#'and \code{M} containing the LMS parameters.}
#'\item{list("\"BCCG\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma} and \code{nu} containing the Box-Cox Cole-Green parameters.}
#'\item{list("\"BCPE\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma}, \code{nu} and \code{tau} containing the Box-Cox Power
#'Exponential parameters.} \item{list("\"BCT\"")}{\code{ref} should contain
#'columns \code{mu}, \code{sigma}, \code{nu} and \code{tau} containing the
#'Box-Cox T distribution parameters.} }
#'
#'@aliases z2y
#'@param z A numerical vector containing standard deviation scores that are to
#'be converted. The length \code{length(z)} determines the size of the output
#'vector.
#'@param x A vector containing the values of the numerical covariate (typically
#'decimal age or height) at which conversion is desired. Values are replicated
#'to match \code{length(y)}.
#'@param sex A character vector indicating whether the male (\code{"M"}) of
#'female (\code{"F"})reference should be used. Values are replicated to match
#'\code{length(y)}.
#'@param sub A character vector indicating the level of the \code{sub} field of
#'the reference standard defined in \code{ref}
#'@param ref A data frame containing a factor \code{sex}, a numerical variable
#'\code{age} containing the tabulated decimal point ages, and two or more
#'numerical variables with reference values. See details.
#'@param dist A string identifying the type of distribution. Values values are:
#'\code{"NO"}, \code{"BCCG"}, \code{"LMS"}, \code{"BCPE"} and \code{"BCT"}.
#'The default is \code{"LMS"}.
#'@param dec A scalar value indicating the number of decimals used to round the
#'value.
#'@param sex.fallback The level of the \code{sex} field used when no match is
#'found. The default \code{sex.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@param sub.fallback The level of the \code{sub} field used when no match is
#'found. The default \code{sub.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@return For \code{y2z()}: A vector with \code{length(y)} elements containing
#'the standard deviation score. For \code{z2y()}: A vector with
#'\code{length(z)} elements containing quantiles.
#'@author Stef van Buuren, 2010
#'@seealso \code{\link{y2z}}
#'@keywords distribution
#'@examples
#'
#'
#'boys <- boys7482
#'
#'# quantile at SD=0 of age 2 years,
#'# height Dutch boys
#'z2y(z=0, x=2)
#'
#'# same for Dutch girls
#'z2y(z=0, x=2, sex="F")
#'
#'# quantile at SD=c(-1,0,1) of age 2 years, BMI Dutch boys
#'z2y(z=c(-1,0,+1), x=2, ref=nl4.bmi)
#'
#'# 0SD line (P50) in kg of weight for age in 5-10 year, Dutch boys
#'z2y(z=rep(0,6), x=5:10, ref=nl4.wgt)
#'
#'# 95th percentile (P95), age 10 years, wfa, Dutch boys
#'z2y(z=qnorm(0.95), x=10, ref=nl4.wgt)
#'
#'# table of P3, P10, P50, P90, P97 of weight for 5-10 year old dutch boys
#'# age per year
#'age <- 5:10
#'p <- c(0.03,0.1,0.5,0.9,0.97)
#'z <- rep(qnorm(p), length(age))
#'x <- rep(age, each=length(p))
#'w <- matrix(z2y(z, x=x, sex="M", ref=nl4.wgt), ncol=length(p),
#' byrow=TRUE)
#'dimnames(w) <- list(age, p)
#'round(w,1)
#'
#'# standard set of Z-scores of weight for all tabulated ages, boys & girls
#'# and three etnicities
#'sds <- c(-2.5, -2, -1, 0, 1, 2, 2.5)
#'age <- nl4.wgt$x
#'z <- rep(sds, times=length(age))
#'x <- rep(age, each=length(sds))
#'sex <- rep(c("M","F"), each=length(z)/2)
#'w <- z2y(z=z, x=x, sex=sex, ref=nl4.wgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(age, sds)
#'data.frame(sub=nl4.wgt$sub,sex=nl4.wgt$sex,round(w,2), row.names=NULL)
#'
#'# P85 of BMI in 5-8 year old Dutch boys and girls
#'e <- expand.grid(age=5:8, sex=c("M","F"))
#'w <- z2y(z=rep(qnorm(0.85),nrow(e)), x=e$age, sex=e$sex, ref=nl4.bmi)
#'w <- matrix(w, nrow=2, byrow=TRUE)
#'dimnames(w) <- list(c("boys","girls"),5:8)
#'w
#'
#'# data transformation of height z-scores to cm-scale
#'z <- c(-1.83, 0.09, 2.33, 0.81, -1.20)
#'x <- c(8.33, 0.23, 19.2, 24.3, 10)
#'sex <- c("M", "M", "F", "M", "F")
#'round(z2y(z=z, x=x, sex=sex, ref=nl4.hgt), 1)
#'
#'# interpolate published height standard
#'# to daily values, days 0-31, boys
#'# on centiles -2SD, 0SD and +2SD
#'days <- 0:31
#'sds <- c(-2, 0, +2)
#'z <- rep(sds, length(days))
#'x <- rep(round(days/365.25,4), each=length(sds))
#'w <- z2y(z, x, sex="M", ref=nl4.hgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(days, sds)
#'w
#'
#'@export
z2y <- function(z = c(-2, 0, 2),
x = 1,
sex = "M",
sub = "N",
ref = get("nl4.hgt"),
dist = "LMS",
dec = 3,
sex.fallback = NA,
sub.fallback = NA)
{
z2y.grp <- function(z, x, ref, dist = "LMS") {
if (dist=="NO") {
check.names(df=ref, needed=c("x","mean","sd"))
mean <- approx(x=ref[,"x"], y=ref[,"mean"], xout=x)$y
sd <- approx(x=ref[,"x"], y=ref[,"sd"], xout=x)$y
return(mean + z*sd)
}
if (dist=="LMS") {
check.names(df=ref, needed=c("x","L","M","S"))
L <- approx(x=ref[,"x"], y=ref[,"L"], xout=x)$y
M <- approx(x=ref[,"x"], y=ref[,"M"], xout=x)$y
S <- approx(x=ref[,"x"], y=ref[,"S"], xout=x)$y
return(ifelse(L>0.01|L<(-0.01),M*(1+L*S*z)^(1/L),M*exp(S*z)))
}
if (dist=="BCCG") {
check.names(df=ref, needed=c("x","nu","mu","sigma"))
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
return(qBCCG(pnorm(z), mu=mu, sigma=sigma, nu=nu))
}
if (dist=="BCPE") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qBCPE(pnorm(z), mu=mu, sigma=sigma, nu=nu, tau=tau))
}
if (dist=="BCT") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qBCT(pnorm(z), mu=mu, sigma=sigma, nu=nu, tau=tau))
}
stop(paste("Reference type", dist, "not implemented."))
}
if (!is.data.frame(ref)) stop("'ref' should be a data frame.")
n <- length(z)
if (n < 1) stop("'z' must have 1 or more values")
if(!is.vector(z)) return(as.numeric(rep(NA, n)))
if(!is.numeric(z)) return(as.numeric(rep(NA, n)))
x <- rep(x, length.out=length(z))
sex <- rep(sex, length.out=length(z))
sub <- rep(sub, length.out=length(z))
dist <- match.arg(dist, choices=c("NO","LMS","BCCG","BCPE","BCT"))
# available levels in ref: sex, sub
lev.sex <- levels(ref$sex[, drop=TRUE])
lev.sub <- levels(ref$sub[, drop=TRUE])
# replace nomatching levels
idx <- is.na(match(sub, lev.sub))
if (any(idx)) {
sub[idx] <- sub.fallback
if (!is.na(sub.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sub.fallback,"'",sep="")
}
idx <- is.na(match(sex, lev.sex))
if (any(idx)) {
sex[idx] <- sex.fallback
if (!is.na(sex.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sex.fallback,"'",sep="")
}
refs <- with(ref,split(ref, f=list(sub, sex)))
xs <- split(x,list(sub, sex))
zs <- split(z,list(sub, sex))
ys <- vector("list",length(zs))
names(ys) <- names(zs)
for(i in 1:length(zs)) {
name <- names(zs)[i]
if(is.null(refs[[name]])) ys[[name]] <- rep(NA,length=length(zs[[name]]))
else ys[[name]] <- z2y.grp(z=zs[[name]], x=xs[[name]], ref=refs[[name]],
dist=dist)
}
y <- unsplit(ys,f=list(sub,sex))
names(y) <- names(z)
return(round(y,dec))
}
#'Converts measurements to standard deviation scores (SDS)
#'
#'Converts measurements into age- and sex-conditional standard deviation score
#'(SDS) using an external reference.
#'
#'Functions \code{z2y()} and \code{y2z()} are the inverse of each other.
#'
#'The argument \code{dist} determines the statistical distribution. The
#'possibilities are as follows: \describe{ \item{list("\"NO\"")}{\code{ref}
#'should contain columns \code{mean} and \code{sd}, containing the mean and the
#'standard deviation in the external reference population.}
#'\item{list("\"LMS\"")}{\code{ref} should contain columns \code{L}, \code{S}
#'and \code{M} containing the LMS parameters.}
#'\item{list("\"BCCG\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma} and \code{nu} containing the Box-Cox Cole-Green parameters.}
#'\item{list("\"BCPE\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma}, \code{nu} and \code{tau} containing the Box-Cox Power
#'Exponential parameters.} \item{list("\"BCT\"")}{\code{ref} should contain
#'columns \code{mu}, \code{sigma}, \code{nu} and \code{tau} containing the
#'Box-Cox T distribution parameters.} }
#'
#'@aliases y2z
#'@param y A numerical vector containing the outcome measurements. The length
#'\code{length(y)} determines the size of the output vector.
#'@param x A vector containing the values of the numerical covariate (typically
#'decimal age or height) at which conversion is desired. Values are replicated
#'to match \code{length(y)}.
#'@param sex A character vector indicating whether the male (\code{"M"}) of
#'female (\code{"F"})reference should be used. Values are replicated to match
#'\code{length(y)}.
#'@param sub A character vector indicating the level of the \code{sub} field of
#'the reference standard defined in \code{ref}
#'@param ref A data frame containing a factor \code{sex}, a numerical variable
#'\code{age} containing the tabulated decimal point ages, and two or more
#'numerical variables with reference values. See details.
#'@param dist A string identifying the type of distribution. Values values are:
#'\code{"NO"}, \code{"BCCG"}, \code{"LMS"}, \code{"BCPE"} and \code{"BCT"}.
#'The default is \code{"LMS"}.
#'@param dec A scalar value indicating the number of decimals used to round the
#'value.
#'@param sex.fallback The level of the \code{sex} field used when no match is
#'found. The default \code{sex.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@param sub.fallback The level of the \code{sub} field used when no match is
#'found. The default \code{sub.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@param tail.adjust Logical. If \code{TRUE} then the WHO method for
#'tail adjustment is applied. The default is \code{FALSE}.
#'@return For \code{y2z()}: A vector with \code{length(y)} elements containing
#'the standard deviation score. For \code{z2y()}: A vector with
#'\code{length(z)} elements containing quantiles.
#'@author Stef van Buuren, 2010
#'@seealso \code{\link{z2y}}
#'@keywords distribution
#'@examples
#'
#'
#'boys <- boys7482
#'
#'# SDS of height 115 cm at age 5 years,
#'# relative to Dutch boys reference
#'y2z(y=115, x=5)
#'
#'# same relative to Dutch girls
#'y2z(y=115, x=5, sex="F")
#'
#'# SDS of IOTF BMI cut-off value for overweight (boys 2-18)
#'# relative to Dutch boys reference
#'cutoff <- c(
#'18.41, 18.15, 17.89, 17.72, 17.55, 17.49, 17.42, 17.49, 17.55, 17.74,
#'17.92, 18.18, 18.44, 18.77, 19.10, 19.47, 19.84, 20.20, 20.55, 20.89,
#'21.22, 21.57, 21.91, 22.27, 22.62, 22.96, 23.29, 23.60, 23.90, 24.18,
#'24.46, 24.73, 25.00)
#'age <- seq(2, 18, by=0.5)
#'(z <- y2z(y=cutoff, x=age, sex="M", ref=nl4.bmi))
#'
#'# apply inverse transformation to check calculations
#'round(z2y(z, age, ref=nl4.bmi), 2)
#'cutoff
#'
#'# calculate percentiles of weight 12 kg at 2 years (boys, girls)
#'100*round(pnorm(y2z(y=c(12,12), x=2, sex=c("M","F"), ref=nl4.wgt)),2)
#'
#'# # percentage of children lighter than 15kg at ages 2-5
#'e <- expand.grid(age=2:5, sex=c("M","F"))
#'z <- y2z(y=rep(15,nrow(e)), x=e$age, sex=e$sex, ref=nl4.wgt)
#'w <- matrix(100*round(pnorm(z),2), nrow=2, byrow=TRUE)
#'dimnames(w) <- list(c("boys","girls"),2:5)
#'w
#'
#'# analysis in Z scale
#'hgt.z <- y2z(y=boys$hgt, x=boys$age, sex="M", ref=nl4.hgt)
#'wgt.z <- y2z(y=boys$wgt, x=boys$age, sex="M", ref=nl4.wgt)
#'plot(hgt.z, wgt.z, col="blue")
#'
#'
#'# z2y
#'
#'# quantile at SD=0 of age 2 years,
#'# height Dutch boys
#'z2y(z=0, x=2)
#'
#'# same for Dutch girls
#'z2y(z=0, x=2, sex="F")
#'
#'# quantile at SD=c(-1,0,1) of age 2 years, BMI Dutch boys
#'z2y(z=c(-1,0,+1), x=2, ref=nl4.bmi)
#'
#'# 0SD line (P50) in kg of weight for age in 5-10 year, Dutch boys
#'z2y(z=rep(0,6), x=5:10, ref=nl4.wgt)
#'
#'# 95th percentile (P95), age 10 years, wfa, Dutch boys
#'z2y(z=qnorm(0.95), x=10, ref=nl4.wgt)
#'
#'# table of P3, P10, P50, P90, P97 of weight for 5-10 year old dutch boys
#'# age per year
#'age <- 5:10
#'p <- c(0.03,0.1,0.5,0.9,0.97)
#'z <- rep(qnorm(p), length(age))
#'x <- rep(age, each=length(p))
#'w <- matrix(z2y(z, x=x, sex="M", ref=nl4.wgt), ncol=length(p),
#' byrow=TRUE)
#'dimnames(w) <- list(age, p)
#'round(w,1)
#'
#'# standard set of Z-scores of weight for all tabulated ages, boys & girls
#'# and three etnicities
#'sds <- c(-2.5, -2, -1, 0, 1, 2, 2.5)
#'age <- nl4.wgt$x
#'z <- rep(sds, times=length(age))
#'x <- rep(age, each=length(sds))
#'sex <- rep(c("M","F"), each=length(z)/2)
#'w <- z2y(z=z, x=x, sex=sex, ref=nl4.wgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(age, sds)
#'data.frame(sub=nl4.wgt$sub,sex=nl4.wgt$sex,round(w,2), row.names=NULL)
#'
#'# P85 of BMI in 5-8 year old Dutch boys and girls
#'e <- expand.grid(age=5:8, sex=c("M","F"))
#'w <- z2y(z=rep(qnorm(0.85),nrow(e)), x=e$age, sex=e$sex, ref=nl4.bmi)
#'w <- matrix(w, nrow=2, byrow=TRUE)
#'dimnames(w) <- list(c("boys","girls"),5:8)
#'w
#'
#'# data transformation of height z-scores to cm-scale
#'z <- c(-1.83, 0.09, 2.33, 0.81, -1.20)
#'x <- c(8.33, 0.23, 19.2, 24.3, 10)
#'sex <- c("M", "M", "F", "M", "F")
#'round(z2y(z=z, x=x, sex=sex, ref=nl4.hgt), 1)
#'
#'# interpolate published height standard
#'# to daily values, days 0-31, boys
#'# on centiles -2SD, 0SD and +2SD
#'days <- 0:31
#'sds <- c(-2, 0, +2)
#'z <- rep(sds, length(days))
#'x <- rep(round(days/365.25,4), each=length(sds))
#'w <- z2y(z, x, sex="M", ref=nl4.hgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(days, sds)
#'w
#'
#'@export
y2z <- function(y = c(75, 80, 85),
x = 1,
sex = "M",
sub = "N",
ref = get("nl4.hgt"),
dist = "LMS",
dec = 3,
sex.fallback = NA,
sub.fallback = NA,
tail.adjust = FALSE)
{
y2z.grp <- function(y, x, ref, dist = "LMS", dec = 3, tail.adjust = FALSE){
if (dist=="NO") {
check.names(df=ref, needed=c("x","mean","sd"))
mean <- approx(x=ref[,"x"], y=ref[,"mean"], xout=x)$y
sd <- approx(x=ref[,"x"], y=ref[,"sd"], xout=x)$y
return((y-mean)/sd)
}
if (dist=="LMS") {
check.names(df=ref, needed=c("x","L","M","S"))
L <- approx(x=ref[,"x"], y=ref[,"L"], xout=x)$y
M <- approx(x=ref[,"x"], y=ref[,"M"], xout=x)$y
S <- approx(x=ref[,"x"], y=ref[,"S"], xout=x)$y
z <- ifelse(L>0.01 | L<(-0.01), (((y/M)^L)-1)/(L*S), log(y/M)/S)
if (tail.adjust) z <- adjust.tail(y, z, L, M, S)
return(z)
}
if (dist=="BCCG") {
check.names(df=ref, needed=c("x","nu","mu","sigma"))
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
return(qnorm(pBCCG(y, mu=mu, sigma=sigma, nu=nu)))
}
if (dist=="BCPE") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qnorm(pBCPE(y, mu=mu, sigma=sigma, nu=nu, tau=tau)))
}
if (dist=="BCT") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qnorm(pBCT(y, mu=mu, sigma=sigma, nu=nu, tau=tau)))
}
if (dist=="BCCG") {
check.names(df=ref, needed=c("x","mu","sigma","nu"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
return()
#return(ifelse(L>0.01|L<(-0.01),(((y/M)^L)-1)/(L*S),log(y/M)/S))
}
stop(paste("Reference type", dist, "not implemented."))
}
if (!is.data.frame(ref)) stop("'ref' should be a data frame.")
n <- length(y)
if (n < 1) stop("'y' must have 1 or more values")
if(!is.vector(y)) return(as.numeric(rep(NA, n)))
if(!is.numeric(y)) return(as.numeric(rep(NA, n)))
x <- rep(x, length.out=length(y))
sex <- rep(sex, length.out=length(y))
sub <- rep(sub, length.out=length(y))
dist <- match.arg(dist, choices=c("NO","LMS","BCCG","BCPE","BCT"))
# available levels in ref: sex, sub
lev.sex <- levels(ref$sex[, drop=TRUE])
lev.sub <- levels(ref$sub[, drop=TRUE])
# replace nomatching levels
idx <- is.na(match(sub, lev.sub))
if (any(idx)) {
sub[idx] <- sub.fallback
if (!is.na(sub.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sub.fallback,"'",sep="")
}
idx <- is.na(match(sex, lev.sex))
if (any(idx)) {
sex[idx] <- sex.fallback
if (!is.na(sex.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sex.fallback,"'",sep="")
}
refs <- with(ref,split(ref, f=list(sub, sex)), drop = TRUE)
xs <- split(x,list(sub, sex), drop = TRUE)
ys <- split(y,list(sub, sex), drop = TRUE)
zs <- vector("list",length(ys))
names(zs) <- names(ys)
if (length(ys) == 0L) return(rep(NA, length(y)))
for(i in 1:length(ys)) {
name <- names(ys)[i]
if(is.null(refs[[name]])) ys[[name]] <- rep(NA,length=length(ys[[name]]))
else zs[[name]] <- y2z.grp(y=ys[[name]], x=xs[[name]], ref=refs[[name]],
dist=dist, tail.adjust = tail.adjust)
}
z <- unsplit(zs,f=list(sub,sex), drop = TRUE)
names(z) <- names(y)
return(round(z, dec))
}
check.names <- function(df, needed){
if (missing(df)) stop("required argument 'df' not found")
if (missing(needed)) stop("required argument 'needed' not found")
notfound <- is.na(match(needed, names(df)))
if (any(notfound)) stop("Not found: ",paste(needed[notfound],collapse=", "))
}
adjust.tail <- function(y, z, L, M, S){
idx <- !is.na(z) & z > 3
if (any(idx)) {
sd3 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*3)^(1/L), M*exp(S*3))
sd2 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*2)^(1/L), M*exp(S*2))
z[idx] <- (3 + (y - sd3)/(sd3 - sd2))[idx]
}
idx <- !is.na(z) & z < (-3)
if (any(idx)) {
sd3 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*(-3))^(1/L), M*exp(S*(-3)))
sd2 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*(-2))^(1/L), M*exp(S*(-2)))
z[idx] <- (-3 + (y - sd3)/(sd2 - sd3))[idx]
}
return(z)
}
Try the AGD package in your browser
Any scripts or data that you put into this service are public.
AGD documentation built on May 29, 2024, 3:19 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions adjust.tail check.names y2z z2y
Documented in y2z z2y
# agd.r
#
# Tools for course Analysis of Growth Data
### SvB 20dec2014
#'Convert standard deviation scores (SDS) to measurements
#'
#'Converts standard deviation score (SDS) into measurements using
#'an age- and sex-conditional external reference.
#'
#'Functions \code{z2y()} and \code{y2z()} are the inverse of each other.
#'
#'The argument \code{dist} determines the statistical distribution. The
#'possibilities are as follows: \describe{ \item{list("\"NO\"")}{\code{ref}
#'should contain columns \code{mean} and \code{sd}, containing the mean and the
#'standard deviation in the external reference population.}
#'\item{list("\"LMS\"")}{\code{ref} should contain columns \code{L}, \code{S}
#'and \code{M} containing the LMS parameters.}
#'\item{list("\"BCCG\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma} and \code{nu} containing the Box-Cox Cole-Green parameters.}
#'\item{list("\"BCPE\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma}, \code{nu} and \code{tau} containing the Box-Cox Power
#'Exponential parameters.} \item{list("\"BCT\"")}{\code{ref} should contain
#'columns \code{mu}, \code{sigma}, \code{nu} and \code{tau} containing the
#'Box-Cox T distribution parameters.} }
#'
#'@aliases z2y
#'@param z A numerical vector containing standard deviation scores that are to
#'be converted. The length \code{length(z)} determines the size of the output
#'vector.
#'@param x A vector containing the values of the numerical covariate (typically
#'decimal age or height) at which conversion is desired. Values are replicated
#'to match \code{length(y)}.
#'@param sex A character vector indicating whether the male (\code{"M"}) of
#'female (\code{"F"})reference should be used. Values are replicated to match
#'\code{length(y)}.
#'@param sub A character vector indicating the level of the \code{sub} field of
#'the reference standard defined in \code{ref}
#'@param ref A data frame containing a factor \code{sex}, a numerical variable
#'\code{age} containing the tabulated decimal point ages, and two or more
#'numerical variables with reference values. See details.
#'@param dist A string identifying the type of distribution. Values values are:
#'\code{"NO"}, \code{"BCCG"}, \code{"LMS"}, \code{"BCPE"} and \code{"BCT"}.
#'The default is \code{"LMS"}.
#'@param dec A scalar value indicating the number of decimals used to round the
#'value.
#'@param sex.fallback The level of the \code{sex} field used when no match is
#'found. The default \code{sex.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@param sub.fallback The level of the \code{sub} field used when no match is
#'found. The default \code{sub.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@return For \code{y2z()}: A vector with \code{length(y)} elements containing
#'the standard deviation score. For \code{z2y()}: A vector with
#'\code{length(z)} elements containing quantiles.
#'@author Stef van Buuren, 2010
#'@seealso \code{\link{y2z}}
#'@keywords distribution
#'@examples
#'
#'
#'boys <- boys7482
#'
#'# quantile at SD=0 of age 2 years,
#'# height Dutch boys
#'z2y(z=0, x=2)
#'
#'# same for Dutch girls
#'z2y(z=0, x=2, sex="F")
#'
#'# quantile at SD=c(-1,0,1) of age 2 years, BMI Dutch boys
#'z2y(z=c(-1,0,+1), x=2, ref=nl4.bmi)
#'
#'# 0SD line (P50) in kg of weight for age in 5-10 year, Dutch boys
#'z2y(z=rep(0,6), x=5:10, ref=nl4.wgt)
#'
#'# 95th percentile (P95), age 10 years, wfa, Dutch boys
#'z2y(z=qnorm(0.95), x=10, ref=nl4.wgt)
#'
#'# table of P3, P10, P50, P90, P97 of weight for 5-10 year old dutch boys
#'# age per year
#'age <- 5:10
#'p <- c(0.03,0.1,0.5,0.9,0.97)
#'z <- rep(qnorm(p), length(age))
#'x <- rep(age, each=length(p))
#'w <- matrix(z2y(z, x=x, sex="M", ref=nl4.wgt), ncol=length(p),
#' byrow=TRUE)
#'dimnames(w) <- list(age, p)
#'round(w,1)
#'
#'# standard set of Z-scores of weight for all tabulated ages, boys & girls
#'# and three etnicities
#'sds <- c(-2.5, -2, -1, 0, 1, 2, 2.5)
#'age <- nl4.wgt$x
#'z <- rep(sds, times=length(age))
#'x <- rep(age, each=length(sds))
#'sex <- rep(c("M","F"), each=length(z)/2)
#'w <- z2y(z=z, x=x, sex=sex, ref=nl4.wgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(age, sds)
#'data.frame(sub=nl4.wgt$sub,sex=nl4.wgt$sex,round(w,2), row.names=NULL)
#'
#'# P85 of BMI in 5-8 year old Dutch boys and girls
#'e <- expand.grid(age=5:8, sex=c("M","F"))
#'w <- z2y(z=rep(qnorm(0.85),nrow(e)), x=e$age, sex=e$sex, ref=nl4.bmi)
#'w <- matrix(w, nrow=2, byrow=TRUE)
#'dimnames(w) <- list(c("boys","girls"),5:8)
#'w
#'
#'# data transformation of height z-scores to cm-scale
#'z <- c(-1.83, 0.09, 2.33, 0.81, -1.20)
#'x <- c(8.33, 0.23, 19.2, 24.3, 10)
#'sex <- c("M", "M", "F", "M", "F")
#'round(z2y(z=z, x=x, sex=sex, ref=nl4.hgt), 1)
#'
#'# interpolate published height standard
#'# to daily values, days 0-31, boys
#'# on centiles -2SD, 0SD and +2SD
#'days <- 0:31
#'sds <- c(-2, 0, +2)
#'z <- rep(sds, length(days))
#'x <- rep(round(days/365.25,4), each=length(sds))
#'w <- z2y(z, x, sex="M", ref=nl4.hgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(days, sds)
#'w
#'
#'@export
z2y <- function(z = c(-2, 0, 2),
x = 1,
sex = "M",
sub = "N",
ref = get("nl4.hgt"),
dist = "LMS",
dec = 3,
sex.fallback = NA,
sub.fallback = NA)
{
z2y.grp <- function(z, x, ref, dist = "LMS") {
if (dist=="NO") {
check.names(df=ref, needed=c("x","mean","sd"))
mean <- approx(x=ref[,"x"], y=ref[,"mean"], xout=x)$y
sd <- approx(x=ref[,"x"], y=ref[,"sd"], xout=x)$y
return(mean + z*sd)
}
if (dist=="LMS") {
check.names(df=ref, needed=c("x","L","M","S"))
L <- approx(x=ref[,"x"], y=ref[,"L"], xout=x)$y
M <- approx(x=ref[,"x"], y=ref[,"M"], xout=x)$y
S <- approx(x=ref[,"x"], y=ref[,"S"], xout=x)$y
return(ifelse(L>0.01|L<(-0.01),M*(1+L*S*z)^(1/L),M*exp(S*z)))
}
if (dist=="BCCG") {
check.names(df=ref, needed=c("x","nu","mu","sigma"))
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
return(qBCCG(pnorm(z), mu=mu, sigma=sigma, nu=nu))
}
if (dist=="BCPE") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qBCPE(pnorm(z), mu=mu, sigma=sigma, nu=nu, tau=tau))
}
if (dist=="BCT") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qBCT(pnorm(z), mu=mu, sigma=sigma, nu=nu, tau=tau))
}
stop(paste("Reference type", dist, "not implemented."))
}
if (!is.data.frame(ref)) stop("'ref' should be a data frame.")
n <- length(z)
if (n < 1) stop("'z' must have 1 or more values")
if(!is.vector(z)) return(as.numeric(rep(NA, n)))
if(!is.numeric(z)) return(as.numeric(rep(NA, n)))
x <- rep(x, length.out=length(z))
sex <- rep(sex, length.out=length(z))
sub <- rep(sub, length.out=length(z))
dist <- match.arg(dist, choices=c("NO","LMS","BCCG","BCPE","BCT"))
# available levels in ref: sex, sub
lev.sex <- levels(ref$sex[, drop=TRUE])
lev.sub <- levels(ref$sub[, drop=TRUE])
# replace nomatching levels
idx <- is.na(match(sub, lev.sub))
if (any(idx)) {
sub[idx] <- sub.fallback
if (!is.na(sub.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sub.fallback,"'",sep="")
}
idx <- is.na(match(sex, lev.sex))
if (any(idx)) {
sex[idx] <- sex.fallback
if (!is.na(sex.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sex.fallback,"'",sep="")
}
refs <- with(ref,split(ref, f=list(sub, sex)))
xs <- split(x,list(sub, sex))
zs <- split(z,list(sub, sex))
ys <- vector("list",length(zs))
names(ys) <- names(zs)
for(i in 1:length(zs)) {
name <- names(zs)[i]
if(is.null(refs[[name]])) ys[[name]] <- rep(NA,length=length(zs[[name]]))
else ys[[name]] <- z2y.grp(z=zs[[name]], x=xs[[name]], ref=refs[[name]],
dist=dist)
}
y <- unsplit(ys,f=list(sub,sex))
names(y) <- names(z)
return(round(y,dec))
}
#'Converts measurements to standard deviation scores (SDS)
#'
#'Converts measurements into age- and sex-conditional standard deviation score
#'(SDS) using an external reference.
#'
#'Functions \code{z2y()} and \code{y2z()} are the inverse of each other.
#'
#'The argument \code{dist} determines the statistical distribution. The
#'possibilities are as follows: \describe{ \item{list("\"NO\"")}{\code{ref}
#'should contain columns \code{mean} and \code{sd}, containing the mean and the
#'standard deviation in the external reference population.}
#'\item{list("\"LMS\"")}{\code{ref} should contain columns \code{L}, \code{S}
#'and \code{M} containing the LMS parameters.}
#'\item{list("\"BCCG\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma} and \code{nu} containing the Box-Cox Cole-Green parameters.}
#'\item{list("\"BCPE\"")}{\code{ref} should contain columns \code{mu},
#'\code{sigma}, \code{nu} and \code{tau} containing the Box-Cox Power
#'Exponential parameters.} \item{list("\"BCT\"")}{\code{ref} should contain
#'columns \code{mu}, \code{sigma}, \code{nu} and \code{tau} containing the
#'Box-Cox T distribution parameters.} }
#'
#'@aliases y2z
#'@param y A numerical vector containing the outcome measurements. The length
#'\code{length(y)} determines the size of the output vector.
#'@param x A vector containing the values of the numerical covariate (typically
#'decimal age or height) at which conversion is desired. Values are replicated
#'to match \code{length(y)}.
#'@param sex A character vector indicating whether the male (\code{"M"}) of
#'female (\code{"F"})reference should be used. Values are replicated to match
#'\code{length(y)}.
#'@param sub A character vector indicating the level of the \code{sub} field of
#'the reference standard defined in \code{ref}
#'@param ref A data frame containing a factor \code{sex}, a numerical variable
#'\code{age} containing the tabulated decimal point ages, and two or more
#'numerical variables with reference values. See details.
#'@param dist A string identifying the type of distribution. Values values are:
#'\code{"NO"}, \code{"BCCG"}, \code{"LMS"}, \code{"BCPE"} and \code{"BCT"}.
#'The default is \code{"LMS"}.
#'@param dec A scalar value indicating the number of decimals used to round the
#'value.
#'@param sex.fallback The level of the \code{sex} field used when no match is
#'found. The default \code{sex.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@param sub.fallback The level of the \code{sub} field used when no match is
#'found. The default \code{sub.fallback=NA} specifies that
#'unmatched entries should receive a \code{NA} value.
#'@param tail.adjust Logical. If \code{TRUE} then the WHO method for
#'tail adjustment is applied. The default is \code{FALSE}.
#'@return For \code{y2z()}: A vector with \code{length(y)} elements containing
#'the standard deviation score. For \code{z2y()}: A vector with
#'\code{length(z)} elements containing quantiles.
#'@author Stef van Buuren, 2010
#'@seealso \code{\link{z2y}}
#'@keywords distribution
#'@examples
#'
#'
#'boys <- boys7482
#'
#'# SDS of height 115 cm at age 5 years,
#'# relative to Dutch boys reference
#'y2z(y=115, x=5)
#'
#'# same relative to Dutch girls
#'y2z(y=115, x=5, sex="F")
#'
#'# SDS of IOTF BMI cut-off value for overweight (boys 2-18)
#'# relative to Dutch boys reference
#'cutoff <- c(
#'18.41, 18.15, 17.89, 17.72, 17.55, 17.49, 17.42, 17.49, 17.55, 17.74,
#'17.92, 18.18, 18.44, 18.77, 19.10, 19.47, 19.84, 20.20, 20.55, 20.89,
#'21.22, 21.57, 21.91, 22.27, 22.62, 22.96, 23.29, 23.60, 23.90, 24.18,
#'24.46, 24.73, 25.00)
#'age <- seq(2, 18, by=0.5)
#'(z <- y2z(y=cutoff, x=age, sex="M", ref=nl4.bmi))
#'
#'# apply inverse transformation to check calculations
#'round(z2y(z, age, ref=nl4.bmi), 2)
#'cutoff
#'
#'# calculate percentiles of weight 12 kg at 2 years (boys, girls)
#'100*round(pnorm(y2z(y=c(12,12), x=2, sex=c("M","F"), ref=nl4.wgt)),2)
#'
#'# # percentage of children lighter than 15kg at ages 2-5
#'e <- expand.grid(age=2:5, sex=c("M","F"))
#'z <- y2z(y=rep(15,nrow(e)), x=e$age, sex=e$sex, ref=nl4.wgt)
#'w <- matrix(100*round(pnorm(z),2), nrow=2, byrow=TRUE)
#'dimnames(w) <- list(c("boys","girls"),2:5)
#'w
#'
#'# analysis in Z scale
#'hgt.z <- y2z(y=boys$hgt, x=boys$age, sex="M", ref=nl4.hgt)
#'wgt.z <- y2z(y=boys$wgt, x=boys$age, sex="M", ref=nl4.wgt)
#'plot(hgt.z, wgt.z, col="blue")
#'
#'
#'# z2y
#'
#'# quantile at SD=0 of age 2 years,
#'# height Dutch boys
#'z2y(z=0, x=2)
#'
#'# same for Dutch girls
#'z2y(z=0, x=2, sex="F")
#'
#'# quantile at SD=c(-1,0,1) of age 2 years, BMI Dutch boys
#'z2y(z=c(-1,0,+1), x=2, ref=nl4.bmi)
#'
#'# 0SD line (P50) in kg of weight for age in 5-10 year, Dutch boys
#'z2y(z=rep(0,6), x=5:10, ref=nl4.wgt)
#'
#'# 95th percentile (P95), age 10 years, wfa, Dutch boys
#'z2y(z=qnorm(0.95), x=10, ref=nl4.wgt)
#'
#'# table of P3, P10, P50, P90, P97 of weight for 5-10 year old dutch boys
#'# age per year
#'age <- 5:10
#'p <- c(0.03,0.1,0.5,0.9,0.97)
#'z <- rep(qnorm(p), length(age))
#'x <- rep(age, each=length(p))
#'w <- matrix(z2y(z, x=x, sex="M", ref=nl4.wgt), ncol=length(p),
#' byrow=TRUE)
#'dimnames(w) <- list(age, p)
#'round(w,1)
#'
#'# standard set of Z-scores of weight for all tabulated ages, boys & girls
#'# and three etnicities
#'sds <- c(-2.5, -2, -1, 0, 1, 2, 2.5)
#'age <- nl4.wgt$x
#'z <- rep(sds, times=length(age))
#'x <- rep(age, each=length(sds))
#'sex <- rep(c("M","F"), each=length(z)/2)
#'w <- z2y(z=z, x=x, sex=sex, ref=nl4.wgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(age, sds)
#'data.frame(sub=nl4.wgt$sub,sex=nl4.wgt$sex,round(w,2), row.names=NULL)
#'
#'# P85 of BMI in 5-8 year old Dutch boys and girls
#'e <- expand.grid(age=5:8, sex=c("M","F"))
#'w <- z2y(z=rep(qnorm(0.85),nrow(e)), x=e$age, sex=e$sex, ref=nl4.bmi)
#'w <- matrix(w, nrow=2, byrow=TRUE)
#'dimnames(w) <- list(c("boys","girls"),5:8)
#'w
#'
#'# data transformation of height z-scores to cm-scale
#'z <- c(-1.83, 0.09, 2.33, 0.81, -1.20)
#'x <- c(8.33, 0.23, 19.2, 24.3, 10)
#'sex <- c("M", "M", "F", "M", "F")
#'round(z2y(z=z, x=x, sex=sex, ref=nl4.hgt), 1)
#'
#'# interpolate published height standard
#'# to daily values, days 0-31, boys
#'# on centiles -2SD, 0SD and +2SD
#'days <- 0:31
#'sds <- c(-2, 0, +2)
#'z <- rep(sds, length(days))
#'x <- rep(round(days/365.25,4), each=length(sds))
#'w <- z2y(z, x, sex="M", ref=nl4.hgt)
#'w <- matrix(w, ncol=length(sds), byrow=TRUE)
#'dimnames(w) <- list(days, sds)
#'w
#'
#'@export
y2z <- function(y = c(75, 80, 85),
x = 1,
sex = "M",
sub = "N",
ref = get("nl4.hgt"),
dist = "LMS",
dec = 3,
sex.fallback = NA,
sub.fallback = NA,
tail.adjust = FALSE)
{
y2z.grp <- function(y, x, ref, dist = "LMS", dec = 3, tail.adjust = FALSE){
if (dist=="NO") {
check.names(df=ref, needed=c("x","mean","sd"))
mean <- approx(x=ref[,"x"], y=ref[,"mean"], xout=x)$y
sd <- approx(x=ref[,"x"], y=ref[,"sd"], xout=x)$y
return((y-mean)/sd)
}
if (dist=="LMS") {
check.names(df=ref, needed=c("x","L","M","S"))
L <- approx(x=ref[,"x"], y=ref[,"L"], xout=x)$y
M <- approx(x=ref[,"x"], y=ref[,"M"], xout=x)$y
S <- approx(x=ref[,"x"], y=ref[,"S"], xout=x)$y
z <- ifelse(L>0.01 | L<(-0.01), (((y/M)^L)-1)/(L*S), log(y/M)/S)
if (tail.adjust) z <- adjust.tail(y, z, L, M, S)
return(z)
}
if (dist=="BCCG") {
check.names(df=ref, needed=c("x","nu","mu","sigma"))
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
return(qnorm(pBCCG(y, mu=mu, sigma=sigma, nu=nu)))
}
if (dist=="BCPE") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qnorm(pBCPE(y, mu=mu, sigma=sigma, nu=nu, tau=tau)))
}
if (dist=="BCT") {
check.names(df=ref, needed=c("x","nu","mu","sigma","tau"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
tau <- approx(x=ref[,"x"], y=ref[,"tau"], xout=x)$y
return(qnorm(pBCT(y, mu=mu, sigma=sigma, nu=nu, tau=tau)))
}
if (dist=="BCCG") {
check.names(df=ref, needed=c("x","mu","sigma","nu"))
mu <- approx(x=ref[,"x"], y=ref[,"mu"], xout=x)$y
sigma <- approx(x=ref[,"x"], y=ref[,"sigma"], xout=x)$y
nu <- approx(x=ref[,"x"], y=ref[,"nu"], xout=x)$y
return()
#return(ifelse(L>0.01|L<(-0.01),(((y/M)^L)-1)/(L*S),log(y/M)/S))
}
stop(paste("Reference type", dist, "not implemented."))
}
if (!is.data.frame(ref)) stop("'ref' should be a data frame.")
n <- length(y)
if (n < 1) stop("'y' must have 1 or more values")
if(!is.vector(y)) return(as.numeric(rep(NA, n)))
if(!is.numeric(y)) return(as.numeric(rep(NA, n)))
x <- rep(x, length.out=length(y))
sex <- rep(sex, length.out=length(y))
sub <- rep(sub, length.out=length(y))
dist <- match.arg(dist, choices=c("NO","LMS","BCCG","BCPE","BCT"))
# available levels in ref: sex, sub
lev.sex <- levels(ref$sex[, drop=TRUE])
lev.sub <- levels(ref$sub[, drop=TRUE])
# replace nomatching levels
idx <- is.na(match(sub, lev.sub))
if (any(idx)) {
sub[idx] <- sub.fallback
if (!is.na(sub.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sub.fallback,"'",sep="")
}
idx <- is.na(match(sex, lev.sex))
if (any(idx)) {
sex[idx] <- sex.fallback
if (!is.na(sex.fallback))
warning("Entries (n=",sum(idx),") replaced by '",sex.fallback,"'",sep="")
}
refs <- with(ref,split(ref, f=list(sub, sex)), drop = TRUE)
xs <- split(x,list(sub, sex), drop = TRUE)
ys <- split(y,list(sub, sex), drop = TRUE)
zs <- vector("list",length(ys))
names(zs) <- names(ys)
if (length(ys) == 0L) return(rep(NA, length(y)))
for(i in 1:length(ys)) {
name <- names(ys)[i]
if(is.null(refs[[name]])) ys[[name]] <- rep(NA,length=length(ys[[name]]))
else zs[[name]] <- y2z.grp(y=ys[[name]], x=xs[[name]], ref=refs[[name]],
dist=dist, tail.adjust = tail.adjust)
}
z <- unsplit(zs,f=list(sub,sex), drop = TRUE)
names(z) <- names(y)
return(round(z, dec))
}
check.names <- function(df, needed){
if (missing(df)) stop("required argument 'df' not found")
if (missing(needed)) stop("required argument 'needed' not found")
notfound <- is.na(match(needed, names(df)))
if (any(notfound)) stop("Not found: ",paste(needed[notfound],collapse=", "))
}
adjust.tail <- function(y, z, L, M, S){
idx <- !is.na(z) & z > 3
if (any(idx)) {
sd3 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*3)^(1/L), M*exp(S*3))
sd2 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*2)^(1/L), M*exp(S*2))
z[idx] <- (3 + (y - sd3)/(sd3 - sd2))[idx]
}
idx <- !is.na(z) & z < (-3)
if (any(idx)) {
sd3 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*(-3))^(1/L), M*exp(S*(-3)))
sd2 <- ifelse(L>0.01|L<(-0.01), M*(1+L*S*(-2))^(1/L), M*exp(S*(-2)))
z[idx] <- (-3 + (y - sd3)/(sd2 - sd3))[idx]
}
return(z)
}
Try the AGD package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.