R/suffDimReduct2.R
In abnormally-distributed/cvreg: Cross Validation and Robust Estimation Utilities
Defines functions
.gsfun
.glmirwls
pSAVE
pSIR
pPIR
grpQreg
grpTheilSen
grpRobGLM
grpGLM
grpOLS
Documented in
grpGLM
grpOLS
grpQreg
grpRobGLM
grpTheilSen
pPIR
pSAVE
pSIR
#' Grouped Ordinary Least Squares
#'
#' @description This is similar to the gOLS function in the sSDR package, but slightly simplified in usage
#' and better handling of ill-conditioned matrices. Note that the covariates in this method must be
#' numeric, and not grouped dummy variable representing a factor.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable (must be numeric, not categorical)
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#'
#'
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpOLS = function(X, Y, idx, ranks = NULL,tol=1e-5, maxiter=100){
X <- as.matrix(X)
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)){rank<-rep(1,ngroups)}else{rank<-ranks}
grps = idx
groups = as.vector(table(idx))
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc=as.matrix(Scale(Y, center = TRUE, scale = FALSE));Xc=scale(X, center = TRUE, scale = FALSE)
Sigma <- solve(cvreg:::ridgepow(cov(Xc), -1))*((n-1)/n);Tau <- cvreg:::ridgepow(Sigma, -1)
eig=eigen(Sigma);Sigma_half=eig$vectors%*%tcrossprod(diag(sqrt(eig$values)),eig$vectors)
inveig=eigen(Tau);Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (cov(Xc,Yc)*((n-1)/n)); V <- Sigma_half %*% U
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];
coef = tcrossprod(pseudoinverse(crossprod(Xg)),Xg)%*%Yc;
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {
for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}
}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(cvreg:::ridgepow(crossprod(V2),-1),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE)
betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,betahat[[g]])}}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1) {obj_diff=1} else {obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups){G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ 0 + ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = lm(form, cbind.data.frame(y = Y, predictors));log_lik = logLik(model)[[1]];AIC=AIC(model);y.pred = fitted(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpOLS")
return(ans)
}
#' Grouped Generalized Linear Models
#'
#' @description This is similar to the \code{\link{grpOLS}} function, but extended to the case of
#' generalized linear models. Note that the covariates in this method ideally must be numeric,
#' and not grouped dummy variable representing a factor.
#'
#' @param X a model matrix
#' @param Y the outcome variable
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#' @param tol convergence tolerance for IRWLS. Deaults to 1e-8.
#' @param maxiter the maximum number of iterations. defaults to 500.
#' @param family one of "gaussian", "Gamma", "binomial", "poisson", "quasibinomial", "quasipoisson", "inverse.gaussian", or "negative.binomial".
#' The family may also provided as an unquoted evaluation of a family function, ie, 'binomial(link="probit")'.
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpGLM = function(X, Y, idx, family=gaussian(), ranks = NULL, tol=1e-8, maxiter=500){
X <- as.matrix(X)
if (is.character(family)){
if (family=="negative.binomial"){
negbinom <- function(theta = NULL, link="log"){
if (is.null(theta)) theta <- 0
eval(expression({MASS::negative.binomial(theta, link)}))
}
}
family <- switch(family,
"binomial" = quasibinomial(),
"poisson" = quasipoisson(),
"gaussian" = gaussian(),
"gamma" = Gamma(),
"Gamma" = Gamma(),
"quasibinomial" = quasibinomial(),
"quasipoisson" = quasipoisson(),
"inverse.gaussian" = inverse.gaussian(),
"negative.binomial" = negbinom()
)
}
fam <- family$family
if (isFALSE(fam %in% c("gaussian","binomial","poisson","Gamma","quasi","inverse.gaussian","negative.binomial"))) {
stop("Must provide a valid GLM family object.")
}
varfun <- family$var
mu.eta <- family$mu.eta
linkinv <- family$linkinv
linkfun <- family$linkfun
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)) {rank <- rep(1, ngroups)} else {rank <- ranks}
grps = idx
groups = as.vector(table(idx))
if(is.factor(Y) || is.character(Y)){
Y <- as.numeric(as.factor(Y))
if (length(unique(Y)) == 2 && min(Y) == 1) Y <- Y-1
}
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc = Y - mean(Y)
Xc=scale(X, center = TRUE, scale = FALSE)
Sigma <- crossprod(Xc)/n
Tau <- pseudoinverse(Sigma)
eig=eigen(Sigma)
Sigma_half=eig$vectors%*%diag(sqrt(eig$values))%*%t(eig$vectors)
inveig=eigen(Tau)
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (cov(Xc, Yc)*((n-1)/n));V <- Sigma_half %*% U
if (fam == "gaussian"){Yc = as.matrix(Scale(Y, center = TRUE, scale = FALSE))} else{Yc = Y}
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];coef = .glmirwls(cbind(Intercept = rep(1, nrow(Xg)), Xg), as.vector(Y), family = family)[-1];
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(cvreg:::ridgepow(crossprod(V2),-1),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE);betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,betahat[[g]])}}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1){obj_diff=1}else{obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups) {G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ ", paste0("SP", 1:ncol(predictors), collapse = "+")))
if (fam=="negative.binomial"){
model = glm(form, cbind.data.frame(y = Y, predictors), family = quasipoisson())
theta <- MASS::theta.md(Y, model$fitted, model$df.residual)
model = glm(form, cbind.data.frame(y = Y, predictors), family = MASS::negative.binomial(theta))
}
else{
model = glm(form, cbind.data.frame(y = Y, predictors), family = family)
}
linear.predictor = model$linear.predictors;y.pred = fitted(model);log_lik = logLik(model);AIC=AIC(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, linear.predictor = linear.predictor, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpGLM")
return(ans)
}
#' Grouped Generalized Linear Models
#'
#' @description This is similar to the \code{\link{grpOLS}} function, but extended to the case of
#' generalized linear models. This supports the Gaussian, Binomial, Poisson, and Gamma likelihoods.
#' Note that the covariates in this method must be numeric, and not grouped dummy variable representing
#' a factor. The algorithm implemented here differs from the grpOLS function, so the results for the
#' Gaussian likelihood may slightly differ when compared to the grpOLS results.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#' @param tol convergence tolerance for IRWLS. Deaults to 1e-8.
#' @param maxiter the maximum number of iterations. defaults to 500.
#' @param family one of "gaussian", "Gamma", "binomial", "poisson", "quasibinomial", "quasipoisson", "inverse.gaussian", or "negative.binomial".
#' The family may also provided as an unquoted evaluation of a family function, ie, 'binomial(link="probit")'.
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpRobGLM = function(X, Y, idx, family=gaussian(), ranks = NULL, tol=1e-8, maxiter=500){
X <- as.matrix(X)
if (is.character(family)){
family <- switch(family,
"binomial" = binomial(),
"poisson" = poisson(),
"gaussian" = gaussian(),
"gamma" = Gamma(),
"Gamma" = Gamma(),
"quasibinomial" = quasibinomial(),
"quasipoisson" = quasipoisson()
)
}
fam <- family$family
if (isTRUE(fam %in% c("quasi","inverse.gaussian","negative.binomial"))){
return(cat(crayon::red("'quasi', 'inverse.gaussian', and 'negative.binomial' are not supported in robgrpGLM.")))
}
if (isFALSE(fam %in% c("gaussian","binomial","poisson","Gamma"))) {
return(cat(crayon::red("Must provide a valid GLM family object.")))
}
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)) {rank <- rep(1, ngroups)} else {rank <- ranks}
grps = idx
groups = as.vector(table(idx))
if(is.factor(Y) || is.character(Y)){
Y <- as.numeric(as.factor(Y))
if (length(unique(Y)) == 2 && min(Y) == 1) Y <- Y-1
}
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc = Y - median(Y)
Xc=sweep(X, 2, colMedians(X), "-")
Sigma <- cov.bacon(Xc,center0="Medians")$cov;
Tau <- pseudoinverse(Sigma)
eig=eigen(Sigma)
Sigma_half=eig$vectors%*%diag(sqrt(eig$values))%*%t(eig$vectors)
inveig=eigen(Tau)
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (crossprod(Xc, Yc)/n); V <- Sigma_half %*% U
if (fam == "gaussian"){Yc = as.matrix(Scale(Y, center = TRUE, scale = FALSE))} else{Yc = Y}
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];coef = cvreg:::.Mqle(cbind(Intercept = rep(1, nrow(Xg)), Xg), as.vector(Y), family = family)$basis[-1];
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(pseudoinverse(crossprod(V2)),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE);betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,betahat[[g]])}}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1){obj_diff=1}else{obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups) {G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = robglm(form, cbind.data.frame(y = Y, predictors), family = family)
linear.predictor = model$linear.predictor;y.pred = model$fitted;log_lik = logLik(model);AIC=AIC(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, linear.predictor = linear.predictor, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpGLM")
return(ans)
}
#' Grouped Theil-Sen Robust Regression
#'
#' @description This is analagous to the \code{\link{grpOLS}} function, except the robust
#' Theil-Sen estimator replaces ordinary least squares.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable (must be numeric, not categorical)
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#'
#'
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpTheilSen = function(X, Y, idx, ranks = NULL,tol=1e-5, maxiter=100){
X <- as.matrix(X)
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)){rank<-rep(1,ngroups)}else{rank<-ranks}
grps = idx
groups = as.vector(table(idx))
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc=as.matrix(Scale2(Y, median, scale = NULL));
Xc=as.matrix(Scale2(X, median, scale = NULL));
Sigma <- crossprod(Xc)/n
Tau <- pseudoinverse(Sigma)
eig=eigen(Sigma);
Sigma_half=eig$vectors%*%tcrossprod(diag(sqrt(eig$values)),eig$vectors)
inveig=eigen(Tau);
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (crossprod(Xc,Yc)/n); V <- Sigma_half %*% U
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];
coef = as.vector(cvreg:::theilsen(Xg,Yc)$coef)[-1]
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {
for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}
}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(pseudoinverse(crossprod(V2)),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE)
betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1){
for (g in 2:ngroups){
B=Matrix::bdiag(B,betahat[[g]])
}
}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1) {obj_diff=1} else {obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups){
G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))
}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ 0 + ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = lm(form, cbind.data.frame(y = Y, predictors));log_lik = logLik(model)[[1]];AIC=AIC(model);y.pred = fitted(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpOLS")
return(ans)
}
#' Grouped Quantile Regression
#'
#' @description This is analagous to the \code{\link{grpOLS}} function, except quantile
#' regression replaces ordinary least squares.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable (must be numeric, not categorical)
#' @param idx group id labels
#' @param q the quantile of interest. defaults to 0.50.
#' @param ranks an indicator for each group whether the covariates of said group are active.
#'
#'
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpQreg = function(X, Y, idx, ranks = NULL, q = 0.5, tol=1e-5, maxiter=100){
cov.q <- function(x, q){
quant <- apply(x,2,function(i)quantile(i,q))
return(list(cov=crossprod(sweep(as.matrix(x),2,quant,"-"))/nrow(x),center=quant))
}
X <- as.matrix(X)
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)){rank<-rep(1,ngroups)}else{rank<-ranks}
grps = idx
groups = as.vector(table(idx))
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc=as.matrix(sweep(Y,2,colMedians(Y)));
Xc=as.matrix(sweep(X,2,colMedians(X)));
Sigma <- crossprod(Xc)/n
Tau <- pseudoinverse(Sigma, -1)
eig=eigen(Sigma);
Sigma_half=eig$vectors%*%tcrossprod(diag(sqrt(eig$values)),eig$vectors)
inveig=eigen(Tau);
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (crossprod(Xc,Yc)/n); V <- Sigma_half %*% U
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];
coef = cvreg:::qreg(Yc ~ ., cbind.data.frame(Yc=Yc,Xg), q=q)$coef[-1]
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {
for (g in 2:ngroups){
B=Matrix::bdiag(B,b_init[[g]])
}
}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(pseudoinverse(crossprod(V2)),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE)
betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1){
for (g in 2:ngroups){
B=Matrix::bdiag(B,betahat[[g]])
}
}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1) {obj_diff=1} else {obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups){
G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))
}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ 0 + ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = lm(form, cbind.data.frame(y = Y, predictors));log_lik = logLik(model)[[1]];AIC=AIC(model);y.pred = fitted(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpOLS")
return(ans)
}
#' Partial Parametric Inverse Regression
#'
#' @param fo model formula
#' @param cat.fo a formula of the form "~categorical" with no left hand side,
#' with the name of the categorical predictor over which partial subspaces will
#' be integrated.
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param df the tuning parameter used for the parametric inverse regression. defaults to 3.
#' @param sm the smoother function to be used. one of cubic regression splines ("cr", the default),
#' basis splines ("bs"), or natural splines ("ns"), or orthogonal polynomials ("poly").
#' @return an sdr object
#' @export
#'
#'
pPIR = function(fo, cat.fo, data, rank="all", df=3, sm = c("cr", "ns", "bs", "poly")){
sm <- match.arg(sm)
X = model.matrix(fo, data)[,-1];Y = model.frame(fo, data)[,1]
W = model.frame(cat.fo, data)
W <- as.numeric(as.factor(W[,1]))
if (rank == "all") rank = ncol(X)
xc=t(t(X)-apply(X,2,mean));signrt=cvreg:::matpower(cov(X),-1/2)
xstand=xc%*%signrt;ystand=(Y-mean(Y))/sd(Y)
N = length(Y); W = as.numeric(as.factor(W)); Nw = as.vector(table(W))
levs = unique(as.numeric(W))
M = lapply(levs, function(w) cvreg:::pir(xstand[which(W == w),], ystand[which(W == w)], df, sm)$sdrmat)
M = sapply(1:length(M), function(i) M[[i]]*(Nw[i]/N), simplify = FALSE)
M = Reduce('+', M)
eig = eigen(M)
evectors = eig$vectors[,1:rank]
evalues = eig$values[1:rank]
basis = signrt%*%evectors
predictors = X %*% basis
colnames(basis) <- paste0("SP", 1:ncol(predictors))
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
rownames(basis) <- colnames(X)
MX = model.matrix(~ 0 + ., as.data.frame(predictors))
betas = as.vector(cvreg:::ridgepow(crossprod(MX), -1) %*% (crossprod(MX, Y)))
y.pred <- as.vector(MX %*% betas)
return(structure(list(basis = basis, predictors = predictors, fitted = y.pred, evalues = eigen(M)$values, y.obs = Y), class = "sdr", model = "pPIR"))
}
#' Partial Sliced Inverse Regression
#'
#' @param fo model formula
#' @param cat.fo a formula of the form "~categorical" with no left hand side,
#' with the name of the categorical predictor over which partial subspaces will
#' be integrated.
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#' @return an sdr object
#' @export
#'
pSIR = function(fo, cat.fo, data, slices=6, rank="all", ytype=c("numeric", "categorical")){
ytype <- match.arg(ytype)
X = model.matrix(fo, data)[,-1]
Y = model.frame(fo, data)[,1]
W = model.frame(cat.fo, data)
W <- as.numeric(as.factor(W[,1]))
if (rank == "all") rank = ncol(X)
xc=t(t(X)-apply(X,2,mean))
signrt=cvreg:::matpower(cov(X),-1/2)
xstand=xc%*%signrt
if (ytype == "numeric") ystand=(Y-mean(Y))/sd(Y)
if (ytype == "categorical") {
ystand <- as.numeric(as.factor(Y))
if (min(ystand) == 0) ystand <- ystand+ 1
slices <- min(slices, length(unique(ystand)))
}
if (is.factor(Y) || is.character(Y)){
ytype <- "categorical"
ystand <- as.numeric(as.factor(ystand))
}
N = length(Y)
W = as.numeric(as.factor(W))
Nw = as.vector(table(W))
levs = unique(as.numeric(W))
M = lapply(levs, function(w) cvreg:::.sir(xstand[which(W == w),], ystand[which(W == w)], slices, ytype)$sdrmat)
M = sapply(1:length(M), function(i) M[[i]]*(Nw[i]/N), simplify = FALSE)
M = Reduce('+', M)
eig = eigen(M)
evectors = eig$vectors[,1:rank]
evalues = eig$values[1:rank]
basis = signrt%*%evectors
predictors = X%*%basis
colnames(basis) <- paste0("SP", 1:ncol(predictors))
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
rownames(basis) <- colnames(X)
MX = model.matrix(~ 0 + ., as.data.frame(predictors))
betas = as.vector(cvreg:::ridgepow(crossprod(MX), -1) %*% (crossprod(MX, Y)))
y.pred <- as.vector(MX %*% betas)
return(structure(list(basis = basis, predictors = predictors, fitted = y.pred, evalues = eigen(M)$values, y.obs = Y), class = "sdr", model = "pSIR"))
}
#' Partial Sliced Average Variance Estimator
#'
#' @param fo model formula
#' @param cat.fo a formula of the form "~categorical" with no left hand side,
#' with the name of the categorical predictor over which partial subspaces will
#' be integrated.
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 3. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#' @return an sdr object
#' @export
#'
pSAVE = function(fo, cat.fo, data, rank="all", slices=3, ytype=c("numeric", "categorical")){
ytype <- match.arg(ytype)
X = model.matrix(fo, data)[,-1]
Y = model.frame(fo, data)[,1]
W = model.frame(cat.fo, data)
W <- as.numeric(as.factor(W[,1]))
if (rank == "all") rank = ncol(X)
xc=t(t(X)-apply(X,2,mean))
signrt=cvreg:::matpower(cov(X),-1/2)
xstand=xc%*%signrt
if (ytype == "numeric") ystand=(Y-mean(Y))/sd(Y)
if (ytype == "categorical") {
ystand <- as.numeric(as.factor(Y))
if (min(ystand) == 0) ystand <- ystand+ 1
slices <- min(slices, length(unique(ystand)))
}
if (is.factor(Y) || is.character(Y)){
ytype <- "categorical"
ystand <- as.numeric(as.factor(ystand))
}
N = length(Y);W = as.numeric(as.factor(W));Nw = as.vector(table(W))
levs = unique(as.numeric(W))
M = lapply(levs, function(w) cvreg:::save(xstand[which(W == w),], ystand[which(W == w)], slices, ytype)$sdrmat)
M = sapply(1:length(M), function(i) M[[i]]*(Nw[i]/N), simplify = FALSE)
M = Reduce('+', M)
eig = eigen(M)
evectors = eig$vectors[,1:rank]
evalues = eig$values[1:rank]
basis = signrt%*%evectors
predictors = X%*%basis
colnames(basis) <- paste0("SP", 1:ncol(predictors))
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
rownames(basis) <- colnames(X)
MX = model.matrix(~ 0 + ., as.data.frame(predictors))
betas = as.vector(cvreg:::ridgepow(crossprod(MX), -1) %*% (crossprod(MX, Y)))
y.pred <- as.vector(MX %*% betas)
return(structure(list(basis = basis, predictors = predictors, fitted = y.pred, evalues = eigen(M)$values, y.obs = Y), class = "sdr", model = "pSAVE"))
}
.glmirwls <- function(Xg, y, family, inits = NULL, threshold = 1e-8){
if (is.character(family)){
if (family=="negative.binomial" && is.null(theta)){
nb <- function() {function(theta) eval(expression({MASS::negative.binomial(theta)}))}
yhat <- .biasreducedglm(Xg,y,family = poisson())$fitted
theta <- MASS::theta.md(y, yhat, nrow(Xg)-ncol(Xg))
family <- nb()
family <- family(theta)
}
else{
family <- switch(family,
"binomial" = quasibinomial(),
"poisson" = quasipoisson(),
"gaussian" = gaussian(),
"gamma" = Gamma(),
"Gamma" = Gamma(),
"quasibinomial" = quasibinomial(),
"quasipoisson" = quasipoisson(),
"inverse.gaussian" = inverse.gaussian(),
"negative.binomial" = eval(expression({MASS::negative.binomial(theta)}))
)
}
}
varfun <- family$var
mu.eta <- family$mu.eta
linkinv <- family$linkinv
linkfun <- family$linkfun
max_iter = 1000;diff = 10000;iter_count = 0
calc_b = function(X,beta){beta = as.vector(beta);return(drop(X%*%beta))}
if (is.null(inits)) {beta=rep(0, ncol(Xg))}else{beta= inits}
while(diff > threshold & iter_count < max_iter) {
eta = as.vector(calc_b(Xg, beta))
p = as.vector(linkinv(eta))
W = diag(mu.eta(p)^2/as.vector(varfun(p)))
xtxi <- pseudoinverse(crossprod(Xg))
beta_change = xtxi%*%crossprod(Xg, y - p)
beta = beta + beta_change
diff = sum(beta_change^2)
iter_count = iter_count + 1
}
return(beta)
}
.gsfun <- function(r, b){
if (r == 0){
Matrix::Matrix(as.matrix(b))
}
else{
Matrix::Matrix(gram.schmidt(as.matrix(b), tol = 1e-100)$Q)
}
}
abnormally-distributed/cvreg documentation built on May 3, 2020, 3:45 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 .gsfun .glmirwls pSAVE pSIR pPIR grpQreg grpTheilSen grpRobGLM grpGLM grpOLS
Documented in grpGLM grpOLS grpQreg grpRobGLM grpTheilSen pPIR pSAVE pSIR
#' Grouped Ordinary Least Squares
#'
#' @description This is similar to the gOLS function in the sSDR package, but slightly simplified in usage
#' and better handling of ill-conditioned matrices. Note that the covariates in this method must be
#' numeric, and not grouped dummy variable representing a factor.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable (must be numeric, not categorical)
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#'
#'
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpOLS = function(X, Y, idx, ranks = NULL,tol=1e-5, maxiter=100){
X <- as.matrix(X)
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)){rank<-rep(1,ngroups)}else{rank<-ranks}
grps = idx
groups = as.vector(table(idx))
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc=as.matrix(Scale(Y, center = TRUE, scale = FALSE));Xc=scale(X, center = TRUE, scale = FALSE)
Sigma <- solve(cvreg:::ridgepow(cov(Xc), -1))*((n-1)/n);Tau <- cvreg:::ridgepow(Sigma, -1)
eig=eigen(Sigma);Sigma_half=eig$vectors%*%tcrossprod(diag(sqrt(eig$values)),eig$vectors)
inveig=eigen(Tau);Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (cov(Xc,Yc)*((n-1)/n)); V <- Sigma_half %*% U
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];
coef = tcrossprod(pseudoinverse(crossprod(Xg)),Xg)%*%Yc;
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {
for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}
}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(cvreg:::ridgepow(crossprod(V2),-1),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE)
betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,betahat[[g]])}}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1) {obj_diff=1} else {obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups){G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ 0 + ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = lm(form, cbind.data.frame(y = Y, predictors));log_lik = logLik(model)[[1]];AIC=AIC(model);y.pred = fitted(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpOLS")
return(ans)
}
#' Grouped Generalized Linear Models
#'
#' @description This is similar to the \code{\link{grpOLS}} function, but extended to the case of
#' generalized linear models. Note that the covariates in this method ideally must be numeric,
#' and not grouped dummy variable representing a factor.
#'
#' @param X a model matrix
#' @param Y the outcome variable
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#' @param tol convergence tolerance for IRWLS. Deaults to 1e-8.
#' @param maxiter the maximum number of iterations. defaults to 500.
#' @param family one of "gaussian", "Gamma", "binomial", "poisson", "quasibinomial", "quasipoisson", "inverse.gaussian", or "negative.binomial".
#' The family may also provided as an unquoted evaluation of a family function, ie, 'binomial(link="probit")'.
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpGLM = function(X, Y, idx, family=gaussian(), ranks = NULL, tol=1e-8, maxiter=500){
X <- as.matrix(X)
if (is.character(family)){
if (family=="negative.binomial"){
negbinom <- function(theta = NULL, link="log"){
if (is.null(theta)) theta <- 0
eval(expression({MASS::negative.binomial(theta, link)}))
}
}
family <- switch(family,
"binomial" = quasibinomial(),
"poisson" = quasipoisson(),
"gaussian" = gaussian(),
"gamma" = Gamma(),
"Gamma" = Gamma(),
"quasibinomial" = quasibinomial(),
"quasipoisson" = quasipoisson(),
"inverse.gaussian" = inverse.gaussian(),
"negative.binomial" = negbinom()
)
}
fam <- family$family
if (isFALSE(fam %in% c("gaussian","binomial","poisson","Gamma","quasi","inverse.gaussian","negative.binomial"))) {
stop("Must provide a valid GLM family object.")
}
varfun <- family$var
mu.eta <- family$mu.eta
linkinv <- family$linkinv
linkfun <- family$linkfun
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)) {rank <- rep(1, ngroups)} else {rank <- ranks}
grps = idx
groups = as.vector(table(idx))
if(is.factor(Y) || is.character(Y)){
Y <- as.numeric(as.factor(Y))
if (length(unique(Y)) == 2 && min(Y) == 1) Y <- Y-1
}
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc = Y - mean(Y)
Xc=scale(X, center = TRUE, scale = FALSE)
Sigma <- crossprod(Xc)/n
Tau <- pseudoinverse(Sigma)
eig=eigen(Sigma)
Sigma_half=eig$vectors%*%diag(sqrt(eig$values))%*%t(eig$vectors)
inveig=eigen(Tau)
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (cov(Xc, Yc)*((n-1)/n));V <- Sigma_half %*% U
if (fam == "gaussian"){Yc = as.matrix(Scale(Y, center = TRUE, scale = FALSE))} else{Yc = Y}
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];coef = .glmirwls(cbind(Intercept = rep(1, nrow(Xg)), Xg), as.vector(Y), family = family)[-1];
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(cvreg:::ridgepow(crossprod(V2),-1),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE);betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,betahat[[g]])}}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1){obj_diff=1}else{obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups) {G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ ", paste0("SP", 1:ncol(predictors), collapse = "+")))
if (fam=="negative.binomial"){
model = glm(form, cbind.data.frame(y = Y, predictors), family = quasipoisson())
theta <- MASS::theta.md(Y, model$fitted, model$df.residual)
model = glm(form, cbind.data.frame(y = Y, predictors), family = MASS::negative.binomial(theta))
}
else{
model = glm(form, cbind.data.frame(y = Y, predictors), family = family)
}
linear.predictor = model$linear.predictors;y.pred = fitted(model);log_lik = logLik(model);AIC=AIC(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, linear.predictor = linear.predictor, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpGLM")
return(ans)
}
#' Grouped Generalized Linear Models
#'
#' @description This is similar to the \code{\link{grpOLS}} function, but extended to the case of
#' generalized linear models. This supports the Gaussian, Binomial, Poisson, and Gamma likelihoods.
#' Note that the covariates in this method must be numeric, and not grouped dummy variable representing
#' a factor. The algorithm implemented here differs from the grpOLS function, so the results for the
#' Gaussian likelihood may slightly differ when compared to the grpOLS results.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#' @param tol convergence tolerance for IRWLS. Deaults to 1e-8.
#' @param maxiter the maximum number of iterations. defaults to 500.
#' @param family one of "gaussian", "Gamma", "binomial", "poisson", "quasibinomial", "quasipoisson", "inverse.gaussian", or "negative.binomial".
#' The family may also provided as an unquoted evaluation of a family function, ie, 'binomial(link="probit")'.
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpRobGLM = function(X, Y, idx, family=gaussian(), ranks = NULL, tol=1e-8, maxiter=500){
X <- as.matrix(X)
if (is.character(family)){
family <- switch(family,
"binomial" = binomial(),
"poisson" = poisson(),
"gaussian" = gaussian(),
"gamma" = Gamma(),
"Gamma" = Gamma(),
"quasibinomial" = quasibinomial(),
"quasipoisson" = quasipoisson()
)
}
fam <- family$family
if (isTRUE(fam %in% c("quasi","inverse.gaussian","negative.binomial"))){
return(cat(crayon::red("'quasi', 'inverse.gaussian', and 'negative.binomial' are not supported in robgrpGLM.")))
}
if (isFALSE(fam %in% c("gaussian","binomial","poisson","Gamma"))) {
return(cat(crayon::red("Must provide a valid GLM family object.")))
}
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)) {rank <- rep(1, ngroups)} else {rank <- ranks}
grps = idx
groups = as.vector(table(idx))
if(is.factor(Y) || is.character(Y)){
Y <- as.numeric(as.factor(Y))
if (length(unique(Y)) == 2 && min(Y) == 1) Y <- Y-1
}
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc = Y - median(Y)
Xc=sweep(X, 2, colMedians(X), "-")
Sigma <- cov.bacon(Xc,center0="Medians")$cov;
Tau <- pseudoinverse(Sigma)
eig=eigen(Sigma)
Sigma_half=eig$vectors%*%diag(sqrt(eig$values))%*%t(eig$vectors)
inveig=eigen(Tau)
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (crossprod(Xc, Yc)/n); V <- Sigma_half %*% U
if (fam == "gaussian"){Yc = as.matrix(Scale(Y, center = TRUE, scale = FALSE))} else{Yc = Y}
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];coef = cvreg:::.Mqle(cbind(Intercept = rep(1, nrow(Xg)), Xg), as.vector(Y), family = family)$basis[-1];
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(pseudoinverse(crossprod(V2)),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE);betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1) {for (g in 2:ngroups) {B=Matrix::bdiag(B,betahat[[g]])}}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1){obj_diff=1}else{obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups) {G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = robglm(form, cbind.data.frame(y = Y, predictors), family = family)
linear.predictor = model$linear.predictor;y.pred = model$fitted;log_lik = logLik(model);AIC=AIC(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, linear.predictor = linear.predictor, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpGLM")
return(ans)
}
#' Grouped Theil-Sen Robust Regression
#'
#' @description This is analagous to the \code{\link{grpOLS}} function, except the robust
#' Theil-Sen estimator replaces ordinary least squares.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable (must be numeric, not categorical)
#' @param idx group id labels
#' @param ranks an indicator for each group whether the covariates of said group are active.
#'
#'
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpTheilSen = function(X, Y, idx, ranks = NULL,tol=1e-5, maxiter=100){
X <- as.matrix(X)
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)){rank<-rep(1,ngroups)}else{rank<-ranks}
grps = idx
groups = as.vector(table(idx))
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc=as.matrix(Scale2(Y, median, scale = NULL));
Xc=as.matrix(Scale2(X, median, scale = NULL));
Sigma <- crossprod(Xc)/n
Tau <- pseudoinverse(Sigma)
eig=eigen(Sigma);
Sigma_half=eig$vectors%*%tcrossprod(diag(sqrt(eig$values)),eig$vectors)
inveig=eigen(Tau);
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (crossprod(Xc,Yc)/n); V <- Sigma_half %*% U
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];
coef = as.vector(cvreg:::theilsen(Xg,Yc)$coef)[-1]
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {
for (g in 2:ngroups) {B=Matrix::bdiag(B,b_init[[g]])}
}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(pseudoinverse(crossprod(V2)),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE)
betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1){
for (g in 2:ngroups){
B=Matrix::bdiag(B,betahat[[g]])
}
}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1) {obj_diff=1} else {obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups){
G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))
}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ 0 + ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = lm(form, cbind.data.frame(y = Y, predictors));log_lik = logLik(model)[[1]];AIC=AIC(model);y.pred = fitted(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpOLS")
return(ans)
}
#' Grouped Quantile Regression
#'
#' @description This is analagous to the \code{\link{grpOLS}} function, except quantile
#' regression replaces ordinary least squares.
#'
#' @param X a model matrix (must be numeric, not categorical)
#' @param Y the outcome variable (must be numeric, not categorical)
#' @param idx group id labels
#' @param q the quantile of interest. defaults to 0.50.
#' @param ranks an indicator for each group whether the covariates of said group are active.
#'
#'
#' @return an sdr object
#' @export
#'
#' @references Liu, Y., Chiaromonte, F. and Li, B. (2017) Structured Ordinary Least Squares: A Sufficient Dimension Reduction approach for regressions with partitioned predictors and heterogeneous units. Biom, 73: 529-539. doi:10.1111/biom.12579 \cr
#'
grpQreg = function(X, Y, idx, ranks = NULL, q = 0.5, tol=1e-5, maxiter=100){
cov.q <- function(x, q){
quant <- apply(x,2,function(i)quantile(i,q))
return(list(cov=crossprod(sweep(as.matrix(x),2,quant,"-"))/nrow(x),center=quant))
}
X <- as.matrix(X)
ngroups = length(unique(idx))
grp.list = lapply(unique(idx), function(i) X[, which(idx == i)])
varnames = unlist(sapply(unique(idx), function(i) colnames(X)[which(idx == i)]))
X = as.matrix(do.call(cbind, grp.list))
colnames(X) <- varnames
if (is.null(ranks)){rank<-rep(1,ngroups)}else{rank<-ranks}
grps = idx
groups = as.vector(table(idx))
Y <- as.matrix(Y)
n=nrow(X)
p=ncol(X)
Gamma = diag(p)
ngroups=length(groups)
ggamma=NULL
for (g in 1:ngroups) {
first=sum(groups[0:(g-1)])+1
last=sum(groups[1:g])
tmp= Gamma[,first:last]
if (is.null(ncol(tmp))) tmp=matrix(tmp,nrow=p)
ggamma=c(ggamma,list(tmp))
}
Yc=as.matrix(sweep(Y,2,colMedians(Y)));
Xc=as.matrix(sweep(X,2,colMedians(X)));
Sigma <- crossprod(Xc)/n
Tau <- pseudoinverse(Sigma, -1)
eig=eigen(Sigma);
Sigma_half=eig$vectors%*%tcrossprod(diag(sqrt(eig$values)),eig$vectors)
inveig=eigen(Tau);
Tau_half=inveig$vectors%*%tcrossprod(diag(sqrt(inveig$values)), inveig$vectors)
U <- Tau %*% (crossprod(Xc,Yc)/n); V <- Sigma_half %*% U
b_init=NULL
for (g in 1:ngroups){
if (rank[g] > 0) {
Xg=Xc%*%ggamma[[g]];
coef = cvreg:::qreg(Yc ~ ., cbind.data.frame(Yc=Yc,Xg), q=q)$coef[-1]
b_init = c(b_init,list(coef))
}
else{
coef = matrix(0, nrow=groups[g]);b_init=c(b_init,list(coef))
}
}
B=b_init[[1]]
if (ngroups>1) {
for (g in 2:ngroups){
B=Matrix::bdiag(B,b_init[[g]])
}
}
obj_diff=100;obj_rec=NULL;iter=0
while (obj_diff>tol && iter<maxiter) {
iter=iter+1
A <- as.matrix(Sigma_half %*% B)
C <- as.matrix(tcrossprod(pseudoinverse(crossprod(A)),A)%*%V)
V2=NULL
for (g in 1:ngroups) {
first=sum(rank[0:(g-1)])+1;last=sum(rank[1:g]);fgs=C[first:last,]
if (is.null(ncol(fgs))) fgs=matrix(fgs,ncol=1)
V2=cbind(V2, kronecker(t(fgs),Sigma_half%*%ggamma[[g]]))
}
vec_b=tcrossprod(pseudoinverse(crossprod(V2)),V2)%*%V
betahat=NULL
for (g in 1:ngroups){
if (rank[g] > 0){
pg=groups[g];first=sum((groups*rank)[0:(g-1)])+1;last=sum((groups*rank)[1:g])
vec_bg=vec_b[first:last];bghat=matrix(vec_bg,nrow=pg,byrow=FALSE)
betahat=c(betahat,list(bghat))
}
else{
bghat = matrix(0, nrow=groups[g]);betahat=c(b_init,list(bghat))
}
}
B=betahat[[1]]
if (ngroups>1){
for (g in 2:ngroups){
B=Matrix::bdiag(B,betahat[[g]])
}
}
objfn=sum((Sigma_half%*%U-Sigma_half%*%Gamma%*%B%*%C)^2);obj_rec=c(obj_rec,objfn)
if (iter==1) {obj_diff=1} else {obj_diff=abs(obj_rec[iter-1]-objfn)}}
G = .gsfun(rank[1], as.matrix(betahat[[1]]))
for (g in 2:ngroups){
G = Matrix::bdiag(G, .gsfun(rank[g], as.matrix(betahat[[g]])))
}
predictors <- Xc %*% as.matrix(G)
colnames(G) <- paste0("GRP", 1:ngroups);rownames(G) <- varnames
colnames(B) <- paste0("GRP", 1:ngroups);rownames(B) <- varnames
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
form = as.formula(paste0("y ~ 0 + ", paste0("SP", 1:ncol(predictors), collapse = "+")))
model = lm(form, cbind.data.frame(y = Y, predictors));log_lik = logLik(model)[[1]];AIC=AIC(model);y.pred = fitted(model)
ans=structure(list(coef=B, basis=G, inits = b_init, predictors=predictors, y.pred=y.pred, y.obs = Y, mf = cbind.data.frame(y.obs = Y, X), log_lik = log_lik, AIC = AIC), class = "sdr", model = "grpOLS")
return(ans)
}
#' Partial Parametric Inverse Regression
#'
#' @param fo model formula
#' @param cat.fo a formula of the form "~categorical" with no left hand side,
#' with the name of the categorical predictor over which partial subspaces will
#' be integrated.
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param df the tuning parameter used for the parametric inverse regression. defaults to 3.
#' @param sm the smoother function to be used. one of cubic regression splines ("cr", the default),
#' basis splines ("bs"), or natural splines ("ns"), or orthogonal polynomials ("poly").
#' @return an sdr object
#' @export
#'
#'
pPIR = function(fo, cat.fo, data, rank="all", df=3, sm = c("cr", "ns", "bs", "poly")){
sm <- match.arg(sm)
X = model.matrix(fo, data)[,-1];Y = model.frame(fo, data)[,1]
W = model.frame(cat.fo, data)
W <- as.numeric(as.factor(W[,1]))
if (rank == "all") rank = ncol(X)
xc=t(t(X)-apply(X,2,mean));signrt=cvreg:::matpower(cov(X),-1/2)
xstand=xc%*%signrt;ystand=(Y-mean(Y))/sd(Y)
N = length(Y); W = as.numeric(as.factor(W)); Nw = as.vector(table(W))
levs = unique(as.numeric(W))
M = lapply(levs, function(w) cvreg:::pir(xstand[which(W == w),], ystand[which(W == w)], df, sm)$sdrmat)
M = sapply(1:length(M), function(i) M[[i]]*(Nw[i]/N), simplify = FALSE)
M = Reduce('+', M)
eig = eigen(M)
evectors = eig$vectors[,1:rank]
evalues = eig$values[1:rank]
basis = signrt%*%evectors
predictors = X %*% basis
colnames(basis) <- paste0("SP", 1:ncol(predictors))
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
rownames(basis) <- colnames(X)
MX = model.matrix(~ 0 + ., as.data.frame(predictors))
betas = as.vector(cvreg:::ridgepow(crossprod(MX), -1) %*% (crossprod(MX, Y)))
y.pred <- as.vector(MX %*% betas)
return(structure(list(basis = basis, predictors = predictors, fitted = y.pred, evalues = eigen(M)$values, y.obs = Y), class = "sdr", model = "pPIR"))
}
#' Partial Sliced Inverse Regression
#'
#' @param fo model formula
#' @param cat.fo a formula of the form "~categorical" with no left hand side,
#' with the name of the categorical predictor over which partial subspaces will
#' be integrated.
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 5. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#' @return an sdr object
#' @export
#'
pSIR = function(fo, cat.fo, data, slices=6, rank="all", ytype=c("numeric", "categorical")){
ytype <- match.arg(ytype)
X = model.matrix(fo, data)[,-1]
Y = model.frame(fo, data)[,1]
W = model.frame(cat.fo, data)
W <- as.numeric(as.factor(W[,1]))
if (rank == "all") rank = ncol(X)
xc=t(t(X)-apply(X,2,mean))
signrt=cvreg:::matpower(cov(X),-1/2)
xstand=xc%*%signrt
if (ytype == "numeric") ystand=(Y-mean(Y))/sd(Y)
if (ytype == "categorical") {
ystand <- as.numeric(as.factor(Y))
if (min(ystand) == 0) ystand <- ystand+ 1
slices <- min(slices, length(unique(ystand)))
}
if (is.factor(Y) || is.character(Y)){
ytype <- "categorical"
ystand <- as.numeric(as.factor(ystand))
}
N = length(Y)
W = as.numeric(as.factor(W))
Nw = as.vector(table(W))
levs = unique(as.numeric(W))
M = lapply(levs, function(w) cvreg:::.sir(xstand[which(W == w),], ystand[which(W == w)], slices, ytype)$sdrmat)
M = sapply(1:length(M), function(i) M[[i]]*(Nw[i]/N), simplify = FALSE)
M = Reduce('+', M)
eig = eigen(M)
evectors = eig$vectors[,1:rank]
evalues = eig$values[1:rank]
basis = signrt%*%evectors
predictors = X%*%basis
colnames(basis) <- paste0("SP", 1:ncol(predictors))
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
rownames(basis) <- colnames(X)
MX = model.matrix(~ 0 + ., as.data.frame(predictors))
betas = as.vector(cvreg:::ridgepow(crossprod(MX), -1) %*% (crossprod(MX, Y)))
y.pred <- as.vector(MX %*% betas)
return(structure(list(basis = basis, predictors = predictors, fitted = y.pred, evalues = eigen(M)$values, y.obs = Y), class = "sdr", model = "pSIR"))
}
#' Partial Sliced Average Variance Estimator
#'
#' @param fo model formula
#' @param cat.fo a formula of the form "~categorical" with no left hand side,
#' with the name of the categorical predictor over which partial subspaces will
#' be integrated.
#' @param data a data frame
#' @param rank the desired number of sufficient predictors to return. the default is "all".
#' @param slices the number of slices into which the response variable should be split. defaults to 3. for categorical response variables the maximum allowed is the number of response levels minus one. if set above this, it is silently adjusted.
#' @param ytype either numeric or categorical
#' @return an sdr object
#' @export
#'
pSAVE = function(fo, cat.fo, data, rank="all", slices=3, ytype=c("numeric", "categorical")){
ytype <- match.arg(ytype)
X = model.matrix(fo, data)[,-1]
Y = model.frame(fo, data)[,1]
W = model.frame(cat.fo, data)
W <- as.numeric(as.factor(W[,1]))
if (rank == "all") rank = ncol(X)
xc=t(t(X)-apply(X,2,mean))
signrt=cvreg:::matpower(cov(X),-1/2)
xstand=xc%*%signrt
if (ytype == "numeric") ystand=(Y-mean(Y))/sd(Y)
if (ytype == "categorical") {
ystand <- as.numeric(as.factor(Y))
if (min(ystand) == 0) ystand <- ystand+ 1
slices <- min(slices, length(unique(ystand)))
}
if (is.factor(Y) || is.character(Y)){
ytype <- "categorical"
ystand <- as.numeric(as.factor(ystand))
}
N = length(Y);W = as.numeric(as.factor(W));Nw = as.vector(table(W))
levs = unique(as.numeric(W))
M = lapply(levs, function(w) cvreg:::save(xstand[which(W == w),], ystand[which(W == w)], slices, ytype)$sdrmat)
M = sapply(1:length(M), function(i) M[[i]]*(Nw[i]/N), simplify = FALSE)
M = Reduce('+', M)
eig = eigen(M)
evectors = eig$vectors[,1:rank]
evalues = eig$values[1:rank]
basis = signrt%*%evectors
predictors = X%*%basis
colnames(basis) <- paste0("SP", 1:ncol(predictors))
colnames(predictors) <- paste0("SP", 1:ncol(predictors))
rownames(basis) <- colnames(X)
MX = model.matrix(~ 0 + ., as.data.frame(predictors))
betas = as.vector(cvreg:::ridgepow(crossprod(MX), -1) %*% (crossprod(MX, Y)))
y.pred <- as.vector(MX %*% betas)
return(structure(list(basis = basis, predictors = predictors, fitted = y.pred, evalues = eigen(M)$values, y.obs = Y), class = "sdr", model = "pSAVE"))
}
.glmirwls <- function(Xg, y, family, inits = NULL, threshold = 1e-8){
if (is.character(family)){
if (family=="negative.binomial" && is.null(theta)){
nb <- function() {function(theta) eval(expression({MASS::negative.binomial(theta)}))}
yhat <- .biasreducedglm(Xg,y,family = poisson())$fitted
theta <- MASS::theta.md(y, yhat, nrow(Xg)-ncol(Xg))
family <- nb()
family <- family(theta)
}
else{
family <- switch(family,
"binomial" = quasibinomial(),
"poisson" = quasipoisson(),
"gaussian" = gaussian(),
"gamma" = Gamma(),
"Gamma" = Gamma(),
"quasibinomial" = quasibinomial(),
"quasipoisson" = quasipoisson(),
"inverse.gaussian" = inverse.gaussian(),
"negative.binomial" = eval(expression({MASS::negative.binomial(theta)}))
)
}
}
varfun <- family$var
mu.eta <- family$mu.eta
linkinv <- family$linkinv
linkfun <- family$linkfun
max_iter = 1000;diff = 10000;iter_count = 0
calc_b = function(X,beta){beta = as.vector(beta);return(drop(X%*%beta))}
if (is.null(inits)) {beta=rep(0, ncol(Xg))}else{beta= inits}
while(diff > threshold & iter_count < max_iter) {
eta = as.vector(calc_b(Xg, beta))
p = as.vector(linkinv(eta))
W = diag(mu.eta(p)^2/as.vector(varfun(p)))
xtxi <- pseudoinverse(crossprod(Xg))
beta_change = xtxi%*%crossprod(Xg, y - p)
beta = beta + beta_change
diff = sum(beta_change^2)
iter_count = iter_count + 1
}
return(beta)
}
.gsfun <- function(r, b){
if (r == 0){
Matrix::Matrix(as.matrix(b))
}
else{
Matrix::Matrix(gram.schmidt(as.matrix(b), tol = 1e-100)$Q)
}
}
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.