Nothing
R/goodnessoffit.R
In globaltest: Testing Groups of Covariates/Features for Association with a Response Variable, with Applications to Gene Set Testing
Defines functions
getK
gtKS
gtLI
reweighZ
btensor
bbase
tpower
ndiff
gtPS
Documented in
bbase
btensor
gtKS
gtLI
gtPS
reweighZ
gtPS <- function(response, null, data, model = c("linear", "logistic", "cox", "poisson", "multinomial"), ...,
covs, bdeg = 3, nint= 10, pord = 2, interact = FALSE, robust = FALSE,
termlabels = FALSE, returnZ = FALSE) {
# store the call
call <- match.call()
# missing data
if (missing(data)) data<-NULL
if (is.matrix(data)) data <- as.data.frame(data)
# evaluate response, which may be one of the colnames of data
response <- eval(call$response, data, parent.frame())
# response model
if (is(response,"lm") || is(response,"glm") || is(response,"coxph")) {
data <- eval(response$call$data)
if (is.null(data)) stop("argument \"data\" in the null model is missing")
else data<-as.data.frame(data)
response<-as.formula(response$call$formula)
}
# settle formula null and vector response if response is a formula
if (missing(null))
if (is(response, "formula"))
null <- response
else
stop("argument \"null\" is missing, with no default")
if (is(response, "formula")) {
name.response <- as.character(eval(response)[[2]])
response <- eval(attr(terms(response, data=data), "variables"), data, environment(response))[[attr(terms(response, data=data), "response")]]
} else {
name.response <- deparse(call$response)
}
# covs default
if (missing(covs)) {
if (is(null,"formula")) {
if (termlabels || (null[[length(null)]]==".")) covs<-attr(terms(null, data=data), "term.labels")
else {
covs<-all.vars(null, functions = FALSE, unique=TRUE)
covs<-covs[!(covs%in%name.response)]
}
}
if (is.data.frame(null) || is.vector(null) || is.matrix(null)) {
null <- data.frame(null)
intercept <- apply(null,2,function(x){all(x==1)})
if (all(intercept == FALSE)) covs <- colnames(null)
else covs <- colnames(null)[-intercept]
}
}
# get the model
if (missing(model)) {
if(is(response, "Surv"))
model <- "cox"
else if ((is.factor(response) && length(levels(response)) <= 2) || is.logical(response))
model <- "logistic"
else if (is.factor(response) && length(levels(response)) > 2)
model <- "multinomial"
}
model <- match.arg(tolower(model), c("linear", "logistic", "cox", "poisson", "multinomial"))
# null design X0
if (is(null, "formula")) X0<-model.matrix(null, data) else X0<-null
m <- ncol(X0)
n <- nrow(X0)
form <- formula("response ~ 0 + X0")
if (model=="cox"){
intercept <- apply(X0,2,function(x){all(x==1)})
X0<- X0[,-intercept, drop=F]
}
# make the hat matrix H
if (model=="linear")
H <- X0 %*% solve(crossprod(X0), t(X0))
else {
null.fit <- switch(model,
logistic = glm(form, family = "binomial", data=data),
poisson = glm(form, family = "poisson", data=data),
cox = coxph(form, method = "breslow", data=data),
multinomial = mlogit(form, data=data)
)
if (model=="logistic" || model=="poisson") {
W <- null.fit$weights
sqrtW <- sqrt(W)
X0Whf <- X0 * matrix(sqrtW, n, m)
X0W <- X0 * matrix(W, n, m)
H <- X0 %*% solve(crossprod(X0Whf), t(X0W))
}
else if (model == "cox") {
expci <- exp(null.fit$linear.predictors)
times <- as.vector(null.fit$y)[1:n]
d <- as.vector(null.fit$y)[(n+1):(2*n)]
dtimes <- unique(times[d == 1])
nd <- length(dtimes)
matrixO <- outer(times, dtimes, "==") * matrix(d, n, nd)
ties <- any(colSums(matrixO) > 1)
atrisk <- outer(times, dtimes, ">=")
hazinc <- as.vector(1 / (expci %*% atrisk))
matrixP <- outer(expci, hazinc, "*") * atrisk
matrixPO <- matrixP %*% t(matrixO)
matrixM <- matrix(d, n, nd, byrow=FALSE) * (!atrisk) - matrixPO %*% (!atrisk)
matrixW <- -crossprod(t(matrixPO))
diag(matrixW) <- diag(matrixW) + rowSums(matrixPO)
X0W <- crossprod(X0, matrixW)
X0WX0 <- X0W %*% X0
X0WX0invX0 <- solve(X0WX0, t(X0))
H <- X0 %*% solve(X0WX0, X0W)
}
else if (model == "multinomial") { #temporary
H <- X0 %*% solve(crossprod(X0), t(X0))
}
}
# smooth terms
if (interact) sterms<-paste(paste("s(",paste(covs,collapse=","),")", sep=""), covs)
else sterms<-paste("s(",covs,")", sep="")
# build datacovs
if (is(null,"formula")) {
if (is.null(data)) datacovs<-cbind(model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
else datacovs<-cbind(data,model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
}
else datacovs<-null[,covs]
datacovs<-as.data.frame(as.matrix(datacovs))
covs<-colnames(datacovs)
p<-ncol(datacovs)
# arguments: list
l <- max(length(bdeg)*is.list(bdeg),length(nint)*is.list(nint),length(pord)*is.list(pord))
if ((is.list(bdeg)|is.list(nint)|is.list(pord)) && l>1) {
lnames<-c(names(nint),names(pord),names(bdeg))[1:l]
lvector<-vector("list",l); names(lvector)<-lnames
if (!is.list(bdeg)) bdeg<-lapply(lvector,function(x){x<-bdeg})
if (!is.list(nint)) nint<-lapply(lvector,function(x){x<-nint})
if (!is.list(pord)) pord<-lapply(lvector,function(x){x<-pord})
# check the list of covs
stopifnot(all(sort(names(bdeg))==sort(names(nint))) && all(sort(names(nint))==sort(names(pord))))
bdeg<-lapply(bdeg,function(x){ if (length(x)==1) rep(x,p) else x})
nint<-lapply(nint,function(x){ if (length(x)==1) rep(x,p) else x})
pord<-lapply(pord,function(x){ if (length(x)==1) rep(x,p) else x})
if (interact) lvector<-lapply(lnames,function(x){prod(unlist(nint[x])+unlist(bdeg[x])) - prod(unlist(pord[x])) }) else
lvector<-lapply(lnames,function(x){sum( unlist(nint[x])+unlist(bdeg[x])-unlist(pord[x]) ) } )
names(lvector)<-lnames
cvector<-c(0,cumsum(lvector))
crange<-lapply(1:l,function(x){(cvector[x]+1):(cvector[x+1])})
names(crange)<-lnames
# built Z matrix
if (interact){
Z<-do.call(cbind, lapply(lapply(lapply(lnames, function(x){
(diag(nrow(H))-H) %*% btensor(datacovs, unlist(bdeg[x]), unlist(nint[x]), unlist(pord[x]))
}), function(z){
z*(1/sqrt(sum(diag(crossprod(z)))))
}),unlist))
colnames(Z)<-1:ncol(Z)
}
else {
Z<-do.call(cbind, lapply(
lapply(lnames, function(x){
do.call(cbind, lapply( lapply( lapply(1:p, function(y){
(diag(nrow(H))-H) %*% bbase(datacovs[,y],unlist(bdeg[x])[y],unlist(nint[x])[y]) %*%
t(ndiff(unlist(nint[x])[y]+unlist(bdeg[x])[y],unlist(pord[x])[y])) %*%
solve(crossprod(t( ndiff(unlist(nint[x])[y]+unlist(bdeg[x])[y],unlist(pord[x])[y]))))
}), function(z){
z*(1/sqrt(sum(diag(crossprod(z)))))
}),unlist))
}), unlist))
colnames(Z)<-unlist(lapply(lnames, function(y){
sapply(1:p,function(x){rep(colnames(datacovs)[x],(unlist(nint[y])[x]+unlist(bdeg[y])[x]-unlist(pord[y])[x])) })
}))
}
# result
if (!robust) res<-gt(response=response, alternative=Z, null=null, data=data, model=model,..., subsets=crange )
else res<-gt(response=response, alternative=Z, null=null, data=data, model=model, ... )
res@structure <- lapply(lnames,function(x){data.frame(s.term=sterms,bdeg=unlist(bdeg[x]),nint=unlist(nint[x]),pord=unlist(pord[x]))})
}
# arguments: vector
else {
if (is.list(bdeg)|is.list(nint)|is.list(pord)) {bdeg<-unlist(bdeg);nint<-unlist(nint);pord<-unlist(pord)}
if (length(bdeg)==1) bdeg<-rep(bdeg,p)
if (length(nint)==1) nint<-rep(nint,p)
if (length(pord)==1) pord<-rep(pord,p)
stopifnot(length(c(bdeg,nint,pord))==3*p)
# build Z
if (interact) {
Z<-btensor(datacovs, bdeg, nint, pord)
colnames(Z)<-1:ncol(Z)
}
else {
Z<-do.call(cbind, lapply( lapply( lapply(1:p, function(x){
(diag(nrow(H))-H) %*% bbase(datacovs[,x],bdeg[x],nint[x]) %*% t(ndiff(nint[x]+bdeg[x],pord[x])) %*% solve(crossprod(t( ndiff(nint[x]+bdeg[x],pord[x]))))
}),
function(y){
y*(1/sqrt(sum(diag(crossprod(y)))))
}),unlist))
colnames(Z)<-unlist(sapply(1:p,function(x){paste(colnames(datacovs)[x],1:(nint[x]+bdeg[x]-pord[x]), sep="") }))
}
# result
res<-gt(response=response, alternative=Z, null=null, data=data, model=model,... )
res@structure <- list(data.frame(s.term=sterms,bdeg,nint,pord))
}
if (!returnZ) return(res)
else return(list(result=res, Z=Z))
}
###################################################
# External functions for gtPS
###################################################
### Extract smooth terms from a gt.object object
setGeneric("sterms", function(object, ...) standardGeneric("sterms"))
setMethod("sterms", "gt.object",
function(object) {
if (length(object@structure)==1) as.data.frame(object@structure)
else object@structure
}
)
### Construct difference matrix of order pord
ndiff <- function(n, pord = 2) {
if (pord == 0) D <-diag(n)
else {
if (pord == 1) { D <- diff(diag(n)) }
else { D <- diff(ndiff(n, pord - 1)) }
}
return(D)
}
### Truncated p-th power function
tpower <- function(x, t, p)
(x - t) ^ p * (x > t)
### Construct B-spline basis
bbase <- function(x, bdeg, nint){
xl <- min(x) - 0.001*(max(x)-min(x))
xr <- max(x) + 0.001*(max(x)-min(x))
dx <- (xr - xl) / nint
knots <- seq(xl - bdeg * dx, xr + bdeg * dx, by = dx)
P <- outer(x, knots, tpower, bdeg)
n <- dim(P)[2]
D <- ndiff(n, bdeg + 1) / (gamma(bdeg + 1) * dx ^ bdeg)
B <- (-1) ^ (bdeg + 1) * P %*% t(D)
B
}
### Construct tensor product of B-spline bases (reparameterized by Kroneker sum of penalties)
btensor <- function(xs, bdeg, nint, pord, returnU=FALSE){
p <- ncol(xs)
stopifnot( p > 1 )
qq <- nint + bdeg
stopifnot( prod(qq) < 5000 )
bigB<-do.call(cbind, lapply(lapply(1:p,function(x){
bbase(xs[,x], bdeg[x], nint[x])
}), unlist))
tpB<-bigB[,( cumsum(qq)[1]-qq[1]+1 ):( cumsum(qq)[1] )]
ksP<-t(ndiff(ncol(tpB),pord[1])) %*% ndiff(ncol(tpB),pord[1])
for (j in 2:p) {
crange<-( cumsum(qq)[j]-qq[j]+1 ):(cumsum(qq)[j])
qj<-length(crange)
tpB<-kronecker(tpB, t(rep(1,qj)) ) * kronecker( t(rep(1,ncol(tpB))), bigB[,crange] )
Pj<-t(ndiff(qj,pord[j])) %*% ndiff(qj,pord[j])
ksP<-kronecker(diag(qj),ksP) + kronecker( Pj ,diag(ncol(ksP)) )
}
P <- reparamZ(Z=tpB, K=ksP)
if (!returnU) return(P)
else {
U <- reparamZ(Z=tpB, K=ksP, returnU=TRUE)$U
return(list(P=P, U=U)) }
}
### Reparameterize Z
reparamZ <-function (Z, pord, K=NULL, tol = 1e-10, returnU=FALSE) {
if (is.null(K)) K = t(ndiff(ncol(Z), pord)) %*% ndiff(ncol(Z), pord)
ek <- eigen(K, symmetric = TRUE)
nullvals <- ek$values < tol
if (any(nullvals)) {
nullK <- ek$vectors[, nullvals]
U <- Z %*% nullK
colnames(U)<-1:ncol(U)
L <- t(sqrt(ek$values[!nullvals]) * t(ek$vectors[, !nullvals]))
P <- Z %*% L %*% solve(t(L) %*% L)
colnames(P)<-1:ncol(P)
}
else {
U <- matrix(0, nrow = nrow(Z), ncol = 0)
P <- Z
}
if (!returnU) return(P)
else return(list(P=P, U=U))
}
### Standardize Z
reweighZ <- function(Z, null.fit){
X0 <- model.matrix(null.fit)
m <- ncol(X0)
n <- nrow(X0)
# make the hat matrix H
if (attr(null.fit,"class")[1]=="lm")
H <- X0 %*% solve(crossprod(X0), t(X0))
else if ( (attr(null.fit,"class")[1]=="glm") && (null.fit$family[1]=="binomial" || null.fit$family[1]=="poisson") ) {
W <- null.fit$weights
sqrtW <- sqrt(W)
X0Whf <- X0 * matrix(sqrtW, n, m)
X0W <- X0 * matrix(W, n, m)
H <- X0 %*% solve(crossprod(X0Whf), t(X0W))
}
else if (attr(null.fit,"class")[1]=="coxph") {
expci <- exp(null.fit$linear.predictors)
times <- as.vector(null.fit$y)[1:n]
d <- as.vector(null.fit$y)[(n+1):(2*n)]
dtimes <- unique(times[d == 1])
nd <- length(dtimes)
matrixO <- outer(times, dtimes, "==") * matrix(d, n, nd)
ties <- any(colSums(matrixO) > 1)
atrisk <- outer(times, dtimes, ">=")
hazinc <- as.vector(1 / (expci %*% atrisk))
matrixP <- outer(expci, hazinc, "*") * atrisk
matrixPO <- matrixP %*% t(matrixO)
matrixM <- matrix(d, n, nd, byrow=FALSE) * (!atrisk) - matrixPO %*% (!atrisk)
matrixW <- -crossprod(t(matrixPO))
diag(matrixW) <- diag(matrixW) + rowSums(matrixPO)
X0W <- crossprod(X0, matrixW)
X0WX0 <- X0W %*% X0
X0WX0invX0 <- solve(X0WX0, t(X0))
H <- X0 %*% solve(X0WX0, X0W)
}
Zadj<-(diag(nrow(H))-H) %*% Z
Zstar <- Zadj*(1/sqrt(sum(diag(crossprod(Zadj)))))
return(Zstar)
}
gtLI <- function(response, null, data, ..., covs, iorder=2,
termlabels = FALSE, standardize = FALSE) {
# store the call
call <- match.call()
# missing data
if (missing(data)) data<-NULL
if (is.matrix(data)) data <- as.data.frame(data)
# evaluate response, which may be one of the colnames of data
response <- eval(call$response, data, parent.frame())
# response model
if (is(response,"lm") | is(response,"glm") | is(response,"coxph")) {
data <- eval(response$call$data)
if (is.null(data)) stop("argument \"data\" in the response model is missing") else data<-as.data.frame(data)
response<-as.formula(response$call$formula)
}
# settle formula null and vector response if response is a formula
if (missing(null))
if (is(response, "formula"))
null <- response
else
stop("argument \"null\" is missing, with no default")
if (is(response, "formula")) {
name.response <- as.character(eval(response)[[2]])
response <- eval(attr(terms(response, data=data), "variables"), data, environment(response))[[attr(terms(response, data=data), "response")]]
} else {
name.response <- deparse(call$response)
}
# covs default
if (missing(covs)) {
if (is(null,"formula")) {
if (termlabels || (null[[length(null)]]==".")) covs<-attr(terms(null, data=data), "term.labels")
else {
covs<-all.vars(null, functions = FALSE, unique=TRUE)
covs<-covs[!(covs%in%name.response)]
}
}
if (is.data.frame(null) || is.vector(null) || is.matrix(null)) {
null <- data.frame(null)
intercept <- apply(null,2,function(x){all(x==1)})
if (all(intercept == FALSE)) covs <- colnames(null)
else covs <- colnames(null)[-intercept]
}
}
covsformula<- as.formula( paste("~ ", paste(covs, collapse= "+")) )
# construct alternative formula
interactions <-as.formula(paste("~(", paste(attr( terms(covsformula) , "term.labels"), collapse= "+") , ")^",iorder) )
# remove terms from alternative that are also in null
dup <- attr(terms(interactions, data=data), "term.labels") %in% attr(terms(covsformula, data=data), "term.labels")
if (all(dup)) stop("all covariates in alternative also in null")
if (any(dup)) interactions <- formula(terms(interactions,data=data)[!dup])
if (is(null,"formula")) {
dup2 <- attr(terms(interactions, data=data), "term.labels") %in% attr(terms(null, data=data), "term.labels")
if (all(dup2)) stop("all covariates in the alternative are also in the null")
if (any(dup2)) interactions <- formula(terms(interactions,data=data)[!dup])
}
# result
res<-gt(response=response, alternative=interactions, null=null, data=data, standardize=standardize, ...)
res@structure<-list(interactions=attr(terms.formula(interactions),"term.labels"))
return(res)
}
gtKS <- function(response, null, data, model = c("linear", "logistic", "cox", "poisson", "multinomial"), ...,
covs, quant = .25, metric = c("euclidean", "pearson"),
kernel=c("uniform", "exponential", "triangular","neighbours","gauss"),
robust = FALSE, scale=TRUE, termlabels = FALSE, returnZ = FALSE) {
# store the call
call <- match.call()
metric <- match.arg(metric)
kernel <- match.arg(kernel)
# missing data
if (missing(data)) data<-NULL
if (is.matrix(data)) data <- as.data.frame(data)
# evaluate response, which may be one of the colnames of data
response <- eval(call$response, data, parent.frame())
# response model
if (is(response,"lm") | is(response,"glm") | is(response,"coxph")) {
data <- eval(response$call$data)
if (is.null(data)) stop("argument \"data\" in the response model is missing") else data<-as.data.frame(data)
response<-as.formula(response$call$formula)
}
# settle half formula null and vector response if response is a formula
if (missing(null))
if (is(response, "formula"))
null <- response
else
stop("argument \"null\" is missing, with no default")
if (is(response, "formula")) {
name.response <- as.character(eval(response)[[2]])
response <- eval(attr(terms(response, data=data), "variables"), data, environment(response))[[attr(terms(response, data=data), "response")]]
} else {
name.response <- deparse(call$response)
}
# covs default
if (missing(covs)) {
if (is(null,"formula")) {
if (termlabels || (null[[length(null)]]==".") ) covs<-attr(terms(null, data=data), "term.labels")
else {
covs<-all.vars(null, functions = FALSE, unique=TRUE)
covs<-covs[!(covs%in%name.response)]
}
}
if (is.data.frame(null) || is.vector(null) || is.matrix(null)) {
null <- data.frame(null)
intercept <- apply(null,2,function(x){all(x==1)})
if (all(intercept == FALSE)) covs <- colnames(null)
else covs <- colnames(null)[-intercept]
}
}
# build the matrix of selected covs
if (is(null,"formula")) {
if (is.null(data)) datacovs<-cbind(model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
else datacovs<-cbind(data,model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
}
else datacovs<-null[,covs]
datacovs<-as.data.frame(as.matrix(datacovs))
covs<-colnames(datacovs)
p<-ncol(datacovs)
# get the model
if (missing(model)) {
if(is(response, "Surv"))
model <- "cox"
else if ((is.factor(response) && length(levels(response)) <= 2) || is.logical(response))
model <- "logistic"
else if (is.factor(response) && length(levels(response)) > 2)
model <- "multinomial"
}
model <- match.arg(tolower(model), c("linear", "logistic", "cox", "poisson", "multinomial"))
# null design X0
if (is(null, "formula")) X0<-model.matrix(null, data) else X0<-null
m <- ncol(X0)
n <- nrow(X0)
form <- formula("response ~ 0 + X0")
# make the hat matrix H
if (model=="linear")
H <- X0 %*% solve(crossprod(X0), t(X0))
else {
null.fit <- switch(model,
logistic = glm(form, family = "binomial", data=data),
poisson = glm(form, family = "poisson", data=data),
cox = coxph(form, method = "breslow", data=data),
multinomial = mlogit(form, data=data)
)
if (model=="logistic" || model=="poisson") {
W <- null.fit$weights
sqrtW <- sqrt(W)
X0Whf <- X0 * matrix(sqrtW, n, m)
X0W <- X0 * matrix(W, n, m)
H <- X0 %*% solve(crossprod(X0Whf), t(X0W))
}
else if (model == "cox") {
expci <- exp(null.fit$linear.predictors)
times <- as.vector(null.fit$y)[1:n]
d <- as.vector(null.fit$y)[(n+1):(2*n)]
dtimes <- unique(times[d == 1])
nd <- length(dtimes)
matrixO <- outer(times, dtimes, "==") * matrix(d, n, nd)
ties <- any(colSums(matrixO) > 1)
atrisk <- outer(times, dtimes, ">=")
hazinc <- as.vector(1 / (expci %*% atrisk))
matrixP <- outer(expci, hazinc, "*") * atrisk
matrixPO <- matrixP %*% t(matrixO)
matrixM <- matrix(d, n, nd, byrow=FALSE) * (!atrisk) - matrixPO %*% (!atrisk)
matrixW <- -crossprod(t(matrixPO))
diag(matrixW) <- diag(matrixW) + rowSums(matrixPO)
X0W <- crossprod(X0, matrixW)
X0WX0 <- X0W %*% X0
X0WX0invX0 <- solve(X0WX0, t(X0))
H <- X0 %*% solve(X0WX0, X0W)
}
else if (model == "multinomial") {
H <- X0 %*% solve(crossprod(X0), t(X0))
}
}
# smooths label
smooths<-paste("s(",paste(covs,collapse=","),")", sep="")
# result
l<-length(quant)
Z<- do.call(cbind,lapply(lapply( lapply(1:l, function(x){
(diag(ncol(H))-H) %*% getK(datacovs, covs, quant[x], metric, kernel, scale)
} ), function(x){ x * (1/sqrt(sum(diag(t(x) %*% x)))) }) ,unlist))
colnames(Z)<-1:(l*n)
lvector<-lapply(1:l,function(x){(((x-1)*n)+1):(x*n)})
names(lvector)<-paste("quant",quant)
if (!robust) res<-gt(response=response, alternative=Z, null=null, data=data, ... ,subsets=lvector)
else res<-gt(response=response, alternative=Z, null=null, data=data, ... )
res@structure <- list(data.frame(smooths,quant, metric, kernel))
if (!returnZ) return(res)
else return(list(result=res, Z=Z))
}
###################################################
# External functions for gtKS
###################################################
getK <- function(data, covs, quant=.25,
metric = c("euclidean", "pearson"),
kernel=c("uniform", "exponential", "triangular","neighbours","gauss"),
scale=TRUE ) {
if (!is.data.frame(data)) stop("data should be a data.frame")
if (missing(covs)) covs <- names(data)
if (length(metric)>1) metric <- "euclidean"
if (!all(covs %in% names(data))) stop("incorrect covariate names")
metric <- match.arg(metric)
kernel <- match.arg(kernel)
# metrics (to add: categorical, mixed)
if ( metric == "euclidean" ) {
if (scale) { DD <- as.matrix(dist(scale(data[,covs]))) }
else { DD <- as.matrix(dist(data[,covs])) }
}
if ( metric == "pearson" ) {
DD <- as.matrix(as.dist(1-abs(cor(t(data[,covs])))))
}
# kernels (to add: Epanechnikov, Biweight, Gauss)
if ( kernel == "uniform" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric(abs(x) <= h) } )
}
if ( kernel == "exponential" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric(exp((-x)/h)) } )
}
if ( kernel == "triangular" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric( (h-abs(x))*(abs(x)<=h) ) } )
}
if ( kernel == "neighbours" ) {
U <- apply(DD, 1, function(x) { as.numeric(x <= quantile(x, quant)) } )
}
if ( kernel == "gauss" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric( dnorm(x/h) ) } )
}
U <- sweep(U, 2, colSums(U), "/")
U
}
Try the globaltest package in your browser
Any scripts or data that you put into this service are public.
globaltest documentation built on Nov. 8, 2020, 8:18 p.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 getK gtKS gtLI reweighZ btensor bbase tpower ndiff gtPS
Documented in bbase btensor gtKS gtLI gtPS reweighZ
gtPS <- function(response, null, data, model = c("linear", "logistic", "cox", "poisson", "multinomial"), ...,
covs, bdeg = 3, nint= 10, pord = 2, interact = FALSE, robust = FALSE,
termlabels = FALSE, returnZ = FALSE) {
# store the call
call <- match.call()
# missing data
if (missing(data)) data<-NULL
if (is.matrix(data)) data <- as.data.frame(data)
# evaluate response, which may be one of the colnames of data
response <- eval(call$response, data, parent.frame())
# response model
if (is(response,"lm") || is(response,"glm") || is(response,"coxph")) {
data <- eval(response$call$data)
if (is.null(data)) stop("argument \"data\" in the null model is missing")
else data<-as.data.frame(data)
response<-as.formula(response$call$formula)
}
# settle formula null and vector response if response is a formula
if (missing(null))
if (is(response, "formula"))
null <- response
else
stop("argument \"null\" is missing, with no default")
if (is(response, "formula")) {
name.response <- as.character(eval(response)[[2]])
response <- eval(attr(terms(response, data=data), "variables"), data, environment(response))[[attr(terms(response, data=data), "response")]]
} else {
name.response <- deparse(call$response)
}
# covs default
if (missing(covs)) {
if (is(null,"formula")) {
if (termlabels || (null[[length(null)]]==".")) covs<-attr(terms(null, data=data), "term.labels")
else {
covs<-all.vars(null, functions = FALSE, unique=TRUE)
covs<-covs[!(covs%in%name.response)]
}
}
if (is.data.frame(null) || is.vector(null) || is.matrix(null)) {
null <- data.frame(null)
intercept <- apply(null,2,function(x){all(x==1)})
if (all(intercept == FALSE)) covs <- colnames(null)
else covs <- colnames(null)[-intercept]
}
}
# get the model
if (missing(model)) {
if(is(response, "Surv"))
model <- "cox"
else if ((is.factor(response) && length(levels(response)) <= 2) || is.logical(response))
model <- "logistic"
else if (is.factor(response) && length(levels(response)) > 2)
model <- "multinomial"
}
model <- match.arg(tolower(model), c("linear", "logistic", "cox", "poisson", "multinomial"))
# null design X0
if (is(null, "formula")) X0<-model.matrix(null, data) else X0<-null
m <- ncol(X0)
n <- nrow(X0)
form <- formula("response ~ 0 + X0")
if (model=="cox"){
intercept <- apply(X0,2,function(x){all(x==1)})
X0<- X0[,-intercept, drop=F]
}
# make the hat matrix H
if (model=="linear")
H <- X0 %*% solve(crossprod(X0), t(X0))
else {
null.fit <- switch(model,
logistic = glm(form, family = "binomial", data=data),
poisson = glm(form, family = "poisson", data=data),
cox = coxph(form, method = "breslow", data=data),
multinomial = mlogit(form, data=data)
)
if (model=="logistic" || model=="poisson") {
W <- null.fit$weights
sqrtW <- sqrt(W)
X0Whf <- X0 * matrix(sqrtW, n, m)
X0W <- X0 * matrix(W, n, m)
H <- X0 %*% solve(crossprod(X0Whf), t(X0W))
}
else if (model == "cox") {
expci <- exp(null.fit$linear.predictors)
times <- as.vector(null.fit$y)[1:n]
d <- as.vector(null.fit$y)[(n+1):(2*n)]
dtimes <- unique(times[d == 1])
nd <- length(dtimes)
matrixO <- outer(times, dtimes, "==") * matrix(d, n, nd)
ties <- any(colSums(matrixO) > 1)
atrisk <- outer(times, dtimes, ">=")
hazinc <- as.vector(1 / (expci %*% atrisk))
matrixP <- outer(expci, hazinc, "*") * atrisk
matrixPO <- matrixP %*% t(matrixO)
matrixM <- matrix(d, n, nd, byrow=FALSE) * (!atrisk) - matrixPO %*% (!atrisk)
matrixW <- -crossprod(t(matrixPO))
diag(matrixW) <- diag(matrixW) + rowSums(matrixPO)
X0W <- crossprod(X0, matrixW)
X0WX0 <- X0W %*% X0
X0WX0invX0 <- solve(X0WX0, t(X0))
H <- X0 %*% solve(X0WX0, X0W)
}
else if (model == "multinomial") { #temporary
H <- X0 %*% solve(crossprod(X0), t(X0))
}
}
# smooth terms
if (interact) sterms<-paste(paste("s(",paste(covs,collapse=","),")", sep=""), covs)
else sterms<-paste("s(",covs,")", sep="")
# build datacovs
if (is(null,"formula")) {
if (is.null(data)) datacovs<-cbind(model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
else datacovs<-cbind(data,model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
}
else datacovs<-null[,covs]
datacovs<-as.data.frame(as.matrix(datacovs))
covs<-colnames(datacovs)
p<-ncol(datacovs)
# arguments: list
l <- max(length(bdeg)*is.list(bdeg),length(nint)*is.list(nint),length(pord)*is.list(pord))
if ((is.list(bdeg)|is.list(nint)|is.list(pord)) && l>1) {
lnames<-c(names(nint),names(pord),names(bdeg))[1:l]
lvector<-vector("list",l); names(lvector)<-lnames
if (!is.list(bdeg)) bdeg<-lapply(lvector,function(x){x<-bdeg})
if (!is.list(nint)) nint<-lapply(lvector,function(x){x<-nint})
if (!is.list(pord)) pord<-lapply(lvector,function(x){x<-pord})
# check the list of covs
stopifnot(all(sort(names(bdeg))==sort(names(nint))) && all(sort(names(nint))==sort(names(pord))))
bdeg<-lapply(bdeg,function(x){ if (length(x)==1) rep(x,p) else x})
nint<-lapply(nint,function(x){ if (length(x)==1) rep(x,p) else x})
pord<-lapply(pord,function(x){ if (length(x)==1) rep(x,p) else x})
if (interact) lvector<-lapply(lnames,function(x){prod(unlist(nint[x])+unlist(bdeg[x])) - prod(unlist(pord[x])) }) else
lvector<-lapply(lnames,function(x){sum( unlist(nint[x])+unlist(bdeg[x])-unlist(pord[x]) ) } )
names(lvector)<-lnames
cvector<-c(0,cumsum(lvector))
crange<-lapply(1:l,function(x){(cvector[x]+1):(cvector[x+1])})
names(crange)<-lnames
# built Z matrix
if (interact){
Z<-do.call(cbind, lapply(lapply(lapply(lnames, function(x){
(diag(nrow(H))-H) %*% btensor(datacovs, unlist(bdeg[x]), unlist(nint[x]), unlist(pord[x]))
}), function(z){
z*(1/sqrt(sum(diag(crossprod(z)))))
}),unlist))
colnames(Z)<-1:ncol(Z)
}
else {
Z<-do.call(cbind, lapply(
lapply(lnames, function(x){
do.call(cbind, lapply( lapply( lapply(1:p, function(y){
(diag(nrow(H))-H) %*% bbase(datacovs[,y],unlist(bdeg[x])[y],unlist(nint[x])[y]) %*%
t(ndiff(unlist(nint[x])[y]+unlist(bdeg[x])[y],unlist(pord[x])[y])) %*%
solve(crossprod(t( ndiff(unlist(nint[x])[y]+unlist(bdeg[x])[y],unlist(pord[x])[y]))))
}), function(z){
z*(1/sqrt(sum(diag(crossprod(z)))))
}),unlist))
}), unlist))
colnames(Z)<-unlist(lapply(lnames, function(y){
sapply(1:p,function(x){rep(colnames(datacovs)[x],(unlist(nint[y])[x]+unlist(bdeg[y])[x]-unlist(pord[y])[x])) })
}))
}
# result
if (!robust) res<-gt(response=response, alternative=Z, null=null, data=data, model=model,..., subsets=crange )
else res<-gt(response=response, alternative=Z, null=null, data=data, model=model, ... )
res@structure <- lapply(lnames,function(x){data.frame(s.term=sterms,bdeg=unlist(bdeg[x]),nint=unlist(nint[x]),pord=unlist(pord[x]))})
}
# arguments: vector
else {
if (is.list(bdeg)|is.list(nint)|is.list(pord)) {bdeg<-unlist(bdeg);nint<-unlist(nint);pord<-unlist(pord)}
if (length(bdeg)==1) bdeg<-rep(bdeg,p)
if (length(nint)==1) nint<-rep(nint,p)
if (length(pord)==1) pord<-rep(pord,p)
stopifnot(length(c(bdeg,nint,pord))==3*p)
# build Z
if (interact) {
Z<-btensor(datacovs, bdeg, nint, pord)
colnames(Z)<-1:ncol(Z)
}
else {
Z<-do.call(cbind, lapply( lapply( lapply(1:p, function(x){
(diag(nrow(H))-H) %*% bbase(datacovs[,x],bdeg[x],nint[x]) %*% t(ndiff(nint[x]+bdeg[x],pord[x])) %*% solve(crossprod(t( ndiff(nint[x]+bdeg[x],pord[x]))))
}),
function(y){
y*(1/sqrt(sum(diag(crossprod(y)))))
}),unlist))
colnames(Z)<-unlist(sapply(1:p,function(x){paste(colnames(datacovs)[x],1:(nint[x]+bdeg[x]-pord[x]), sep="") }))
}
# result
res<-gt(response=response, alternative=Z, null=null, data=data, model=model,... )
res@structure <- list(data.frame(s.term=sterms,bdeg,nint,pord))
}
if (!returnZ) return(res)
else return(list(result=res, Z=Z))
}
###################################################
# External functions for gtPS
###################################################
### Extract smooth terms from a gt.object object
setGeneric("sterms", function(object, ...) standardGeneric("sterms"))
setMethod("sterms", "gt.object",
function(object) {
if (length(object@structure)==1) as.data.frame(object@structure)
else object@structure
}
)
### Construct difference matrix of order pord
ndiff <- function(n, pord = 2) {
if (pord == 0) D <-diag(n)
else {
if (pord == 1) { D <- diff(diag(n)) }
else { D <- diff(ndiff(n, pord - 1)) }
}
return(D)
}
### Truncated p-th power function
tpower <- function(x, t, p)
(x - t) ^ p * (x > t)
### Construct B-spline basis
bbase <- function(x, bdeg, nint){
xl <- min(x) - 0.001*(max(x)-min(x))
xr <- max(x) + 0.001*(max(x)-min(x))
dx <- (xr - xl) / nint
knots <- seq(xl - bdeg * dx, xr + bdeg * dx, by = dx)
P <- outer(x, knots, tpower, bdeg)
n <- dim(P)[2]
D <- ndiff(n, bdeg + 1) / (gamma(bdeg + 1) * dx ^ bdeg)
B <- (-1) ^ (bdeg + 1) * P %*% t(D)
B
}
### Construct tensor product of B-spline bases (reparameterized by Kroneker sum of penalties)
btensor <- function(xs, bdeg, nint, pord, returnU=FALSE){
p <- ncol(xs)
stopifnot( p > 1 )
qq <- nint + bdeg
stopifnot( prod(qq) < 5000 )
bigB<-do.call(cbind, lapply(lapply(1:p,function(x){
bbase(xs[,x], bdeg[x], nint[x])
}), unlist))
tpB<-bigB[,( cumsum(qq)[1]-qq[1]+1 ):( cumsum(qq)[1] )]
ksP<-t(ndiff(ncol(tpB),pord[1])) %*% ndiff(ncol(tpB),pord[1])
for (j in 2:p) {
crange<-( cumsum(qq)[j]-qq[j]+1 ):(cumsum(qq)[j])
qj<-length(crange)
tpB<-kronecker(tpB, t(rep(1,qj)) ) * kronecker( t(rep(1,ncol(tpB))), bigB[,crange] )
Pj<-t(ndiff(qj,pord[j])) %*% ndiff(qj,pord[j])
ksP<-kronecker(diag(qj),ksP) + kronecker( Pj ,diag(ncol(ksP)) )
}
P <- reparamZ(Z=tpB, K=ksP)
if (!returnU) return(P)
else {
U <- reparamZ(Z=tpB, K=ksP, returnU=TRUE)$U
return(list(P=P, U=U)) }
}
### Reparameterize Z
reparamZ <-function (Z, pord, K=NULL, tol = 1e-10, returnU=FALSE) {
if (is.null(K)) K = t(ndiff(ncol(Z), pord)) %*% ndiff(ncol(Z), pord)
ek <- eigen(K, symmetric = TRUE)
nullvals <- ek$values < tol
if (any(nullvals)) {
nullK <- ek$vectors[, nullvals]
U <- Z %*% nullK
colnames(U)<-1:ncol(U)
L <- t(sqrt(ek$values[!nullvals]) * t(ek$vectors[, !nullvals]))
P <- Z %*% L %*% solve(t(L) %*% L)
colnames(P)<-1:ncol(P)
}
else {
U <- matrix(0, nrow = nrow(Z), ncol = 0)
P <- Z
}
if (!returnU) return(P)
else return(list(P=P, U=U))
}
### Standardize Z
reweighZ <- function(Z, null.fit){
X0 <- model.matrix(null.fit)
m <- ncol(X0)
n <- nrow(X0)
# make the hat matrix H
if (attr(null.fit,"class")[1]=="lm")
H <- X0 %*% solve(crossprod(X0), t(X0))
else if ( (attr(null.fit,"class")[1]=="glm") && (null.fit$family[1]=="binomial" || null.fit$family[1]=="poisson") ) {
W <- null.fit$weights
sqrtW <- sqrt(W)
X0Whf <- X0 * matrix(sqrtW, n, m)
X0W <- X0 * matrix(W, n, m)
H <- X0 %*% solve(crossprod(X0Whf), t(X0W))
}
else if (attr(null.fit,"class")[1]=="coxph") {
expci <- exp(null.fit$linear.predictors)
times <- as.vector(null.fit$y)[1:n]
d <- as.vector(null.fit$y)[(n+1):(2*n)]
dtimes <- unique(times[d == 1])
nd <- length(dtimes)
matrixO <- outer(times, dtimes, "==") * matrix(d, n, nd)
ties <- any(colSums(matrixO) > 1)
atrisk <- outer(times, dtimes, ">=")
hazinc <- as.vector(1 / (expci %*% atrisk))
matrixP <- outer(expci, hazinc, "*") * atrisk
matrixPO <- matrixP %*% t(matrixO)
matrixM <- matrix(d, n, nd, byrow=FALSE) * (!atrisk) - matrixPO %*% (!atrisk)
matrixW <- -crossprod(t(matrixPO))
diag(matrixW) <- diag(matrixW) + rowSums(matrixPO)
X0W <- crossprod(X0, matrixW)
X0WX0 <- X0W %*% X0
X0WX0invX0 <- solve(X0WX0, t(X0))
H <- X0 %*% solve(X0WX0, X0W)
}
Zadj<-(diag(nrow(H))-H) %*% Z
Zstar <- Zadj*(1/sqrt(sum(diag(crossprod(Zadj)))))
return(Zstar)
}
gtLI <- function(response, null, data, ..., covs, iorder=2,
termlabels = FALSE, standardize = FALSE) {
# store the call
call <- match.call()
# missing data
if (missing(data)) data<-NULL
if (is.matrix(data)) data <- as.data.frame(data)
# evaluate response, which may be one of the colnames of data
response <- eval(call$response, data, parent.frame())
# response model
if (is(response,"lm") | is(response,"glm") | is(response,"coxph")) {
data <- eval(response$call$data)
if (is.null(data)) stop("argument \"data\" in the response model is missing") else data<-as.data.frame(data)
response<-as.formula(response$call$formula)
}
# settle formula null and vector response if response is a formula
if (missing(null))
if (is(response, "formula"))
null <- response
else
stop("argument \"null\" is missing, with no default")
if (is(response, "formula")) {
name.response <- as.character(eval(response)[[2]])
response <- eval(attr(terms(response, data=data), "variables"), data, environment(response))[[attr(terms(response, data=data), "response")]]
} else {
name.response <- deparse(call$response)
}
# covs default
if (missing(covs)) {
if (is(null,"formula")) {
if (termlabels || (null[[length(null)]]==".")) covs<-attr(terms(null, data=data), "term.labels")
else {
covs<-all.vars(null, functions = FALSE, unique=TRUE)
covs<-covs[!(covs%in%name.response)]
}
}
if (is.data.frame(null) || is.vector(null) || is.matrix(null)) {
null <- data.frame(null)
intercept <- apply(null,2,function(x){all(x==1)})
if (all(intercept == FALSE)) covs <- colnames(null)
else covs <- colnames(null)[-intercept]
}
}
covsformula<- as.formula( paste("~ ", paste(covs, collapse= "+")) )
# construct alternative formula
interactions <-as.formula(paste("~(", paste(attr( terms(covsformula) , "term.labels"), collapse= "+") , ")^",iorder) )
# remove terms from alternative that are also in null
dup <- attr(terms(interactions, data=data), "term.labels") %in% attr(terms(covsformula, data=data), "term.labels")
if (all(dup)) stop("all covariates in alternative also in null")
if (any(dup)) interactions <- formula(terms(interactions,data=data)[!dup])
if (is(null,"formula")) {
dup2 <- attr(terms(interactions, data=data), "term.labels") %in% attr(terms(null, data=data), "term.labels")
if (all(dup2)) stop("all covariates in the alternative are also in the null")
if (any(dup2)) interactions <- formula(terms(interactions,data=data)[!dup])
}
# result
res<-gt(response=response, alternative=interactions, null=null, data=data, standardize=standardize, ...)
res@structure<-list(interactions=attr(terms.formula(interactions),"term.labels"))
return(res)
}
gtKS <- function(response, null, data, model = c("linear", "logistic", "cox", "poisson", "multinomial"), ...,
covs, quant = .25, metric = c("euclidean", "pearson"),
kernel=c("uniform", "exponential", "triangular","neighbours","gauss"),
robust = FALSE, scale=TRUE, termlabels = FALSE, returnZ = FALSE) {
# store the call
call <- match.call()
metric <- match.arg(metric)
kernel <- match.arg(kernel)
# missing data
if (missing(data)) data<-NULL
if (is.matrix(data)) data <- as.data.frame(data)
# evaluate response, which may be one of the colnames of data
response <- eval(call$response, data, parent.frame())
# response model
if (is(response,"lm") | is(response,"glm") | is(response,"coxph")) {
data <- eval(response$call$data)
if (is.null(data)) stop("argument \"data\" in the response model is missing") else data<-as.data.frame(data)
response<-as.formula(response$call$formula)
}
# settle half formula null and vector response if response is a formula
if (missing(null))
if (is(response, "formula"))
null <- response
else
stop("argument \"null\" is missing, with no default")
if (is(response, "formula")) {
name.response <- as.character(eval(response)[[2]])
response <- eval(attr(terms(response, data=data), "variables"), data, environment(response))[[attr(terms(response, data=data), "response")]]
} else {
name.response <- deparse(call$response)
}
# covs default
if (missing(covs)) {
if (is(null,"formula")) {
if (termlabels || (null[[length(null)]]==".") ) covs<-attr(terms(null, data=data), "term.labels")
else {
covs<-all.vars(null, functions = FALSE, unique=TRUE)
covs<-covs[!(covs%in%name.response)]
}
}
if (is.data.frame(null) || is.vector(null) || is.matrix(null)) {
null <- data.frame(null)
intercept <- apply(null,2,function(x){all(x==1)})
if (all(intercept == FALSE)) covs <- colnames(null)
else covs <- colnames(null)[-intercept]
}
}
# build the matrix of selected covs
if (is(null,"formula")) {
if (is.null(data)) datacovs<-cbind(model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
else datacovs<-cbind(data,model.frame(null,data=data),model.matrix(null,data=data))[,covs,drop=F]
}
else datacovs<-null[,covs]
datacovs<-as.data.frame(as.matrix(datacovs))
covs<-colnames(datacovs)
p<-ncol(datacovs)
# get the model
if (missing(model)) {
if(is(response, "Surv"))
model <- "cox"
else if ((is.factor(response) && length(levels(response)) <= 2) || is.logical(response))
model <- "logistic"
else if (is.factor(response) && length(levels(response)) > 2)
model <- "multinomial"
}
model <- match.arg(tolower(model), c("linear", "logistic", "cox", "poisson", "multinomial"))
# null design X0
if (is(null, "formula")) X0<-model.matrix(null, data) else X0<-null
m <- ncol(X0)
n <- nrow(X0)
form <- formula("response ~ 0 + X0")
# make the hat matrix H
if (model=="linear")
H <- X0 %*% solve(crossprod(X0), t(X0))
else {
null.fit <- switch(model,
logistic = glm(form, family = "binomial", data=data),
poisson = glm(form, family = "poisson", data=data),
cox = coxph(form, method = "breslow", data=data),
multinomial = mlogit(form, data=data)
)
if (model=="logistic" || model=="poisson") {
W <- null.fit$weights
sqrtW <- sqrt(W)
X0Whf <- X0 * matrix(sqrtW, n, m)
X0W <- X0 * matrix(W, n, m)
H <- X0 %*% solve(crossprod(X0Whf), t(X0W))
}
else if (model == "cox") {
expci <- exp(null.fit$linear.predictors)
times <- as.vector(null.fit$y)[1:n]
d <- as.vector(null.fit$y)[(n+1):(2*n)]
dtimes <- unique(times[d == 1])
nd <- length(dtimes)
matrixO <- outer(times, dtimes, "==") * matrix(d, n, nd)
ties <- any(colSums(matrixO) > 1)
atrisk <- outer(times, dtimes, ">=")
hazinc <- as.vector(1 / (expci %*% atrisk))
matrixP <- outer(expci, hazinc, "*") * atrisk
matrixPO <- matrixP %*% t(matrixO)
matrixM <- matrix(d, n, nd, byrow=FALSE) * (!atrisk) - matrixPO %*% (!atrisk)
matrixW <- -crossprod(t(matrixPO))
diag(matrixW) <- diag(matrixW) + rowSums(matrixPO)
X0W <- crossprod(X0, matrixW)
X0WX0 <- X0W %*% X0
X0WX0invX0 <- solve(X0WX0, t(X0))
H <- X0 %*% solve(X0WX0, X0W)
}
else if (model == "multinomial") {
H <- X0 %*% solve(crossprod(X0), t(X0))
}
}
# smooths label
smooths<-paste("s(",paste(covs,collapse=","),")", sep="")
# result
l<-length(quant)
Z<- do.call(cbind,lapply(lapply( lapply(1:l, function(x){
(diag(ncol(H))-H) %*% getK(datacovs, covs, quant[x], metric, kernel, scale)
} ), function(x){ x * (1/sqrt(sum(diag(t(x) %*% x)))) }) ,unlist))
colnames(Z)<-1:(l*n)
lvector<-lapply(1:l,function(x){(((x-1)*n)+1):(x*n)})
names(lvector)<-paste("quant",quant)
if (!robust) res<-gt(response=response, alternative=Z, null=null, data=data, ... ,subsets=lvector)
else res<-gt(response=response, alternative=Z, null=null, data=data, ... )
res@structure <- list(data.frame(smooths,quant, metric, kernel))
if (!returnZ) return(res)
else return(list(result=res, Z=Z))
}
###################################################
# External functions for gtKS
###################################################
getK <- function(data, covs, quant=.25,
metric = c("euclidean", "pearson"),
kernel=c("uniform", "exponential", "triangular","neighbours","gauss"),
scale=TRUE ) {
if (!is.data.frame(data)) stop("data should be a data.frame")
if (missing(covs)) covs <- names(data)
if (length(metric)>1) metric <- "euclidean"
if (!all(covs %in% names(data))) stop("incorrect covariate names")
metric <- match.arg(metric)
kernel <- match.arg(kernel)
# metrics (to add: categorical, mixed)
if ( metric == "euclidean" ) {
if (scale) { DD <- as.matrix(dist(scale(data[,covs]))) }
else { DD <- as.matrix(dist(data[,covs])) }
}
if ( metric == "pearson" ) {
DD <- as.matrix(as.dist(1-abs(cor(t(data[,covs])))))
}
# kernels (to add: Epanechnikov, Biweight, Gauss)
if ( kernel == "uniform" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric(abs(x) <= h) } )
}
if ( kernel == "exponential" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric(exp((-x)/h)) } )
}
if ( kernel == "triangular" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric( (h-abs(x))*(abs(x)<=h) ) } )
}
if ( kernel == "neighbours" ) {
U <- apply(DD, 1, function(x) { as.numeric(x <= quantile(x, quant)) } )
}
if ( kernel == "gauss" ) {
dDD <- DD[diag(nrow(DD)) == 0]
h <- quantile(dDD, quant)
U <- apply(DD, c(1,2), function(x) { as.numeric( dnorm(x/h) ) } )
}
U <- sweep(U, 2, colSums(U), "/")
U
}
Try the globaltest 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.