Movatterモバイル変換

The Rotated H-transform of a 3d Array by a Matrix

Description

This function is an implementation of the\rho-operator found inCurrie et al 2006. It forms the basis of the GLAM arithmetic.

Usage

RH(M, A)

Arguments

Details

For details seeCurrie et al 2006. Note that this particular implementationis not used in the routines underlying the optimization procedure.

Value

A 3d array of sizep_2 \times p_3 \times n.

Author(s)

Adam Lund

References

Currie, I. D., M. Durban, and P. H. C. Eilers (2006). Generalized lineararray models with applications to multidimensional smoothing.Journal of the Royal Statistical Society. Series B. 68, 259-280. url = http://dx.doi.org/10.1111/j.1467-9868.2006.00543.x.

Examples

n1 <- 65; n2 <- 26; n3 <- 13; p1 <- 13; p2 <- 5; p3 <- 4##marginal design matrices (Kronecker components)X1 <- matrix(rnorm(n1 * p1), n1, p1)X2 <- matrix(rnorm(n2 * p2), n2, p2)X3 <- matrix(rnorm(n3 * p3), n3, p3)Beta <- array(rnorm(p1 * p2 * p3, 0, 1), c(p1 , p2, p3))max(abs(c(RH(X3, RH(X2, RH(X1, Beta)))) - kronecker(X3, kronecker(X2, X1)) %*% c(Beta)))

Soft Maximin Estimation for Large Scale Array Data with Known Groups

Description

Efficient design matrix free procedure for solving a soft maximin problem forlarge scale array-tensor structured models, seeLund et al., 2020.Currently Lasso and SCAD penalized estimation is implemented.

Usage

softmaximin(X,            Y,            zeta,            penalty = c("lasso", "scad"),            alg = c("npg", "fista"),            nlambda = 30,            lambda.min.ratio = 1e-04,            lambda = NULL,            penalty.factor = NULL,            reltol = 1e-05,            maxiter = 15000,            steps = 1,            btmax = 100,            c = 0.0001,            tau = 2,            M = 4,            nu = 1,            Lmin = 0,            log = TRUE)

Arguments

Details

FollowingLund et al., 2020 this package solves the optimization problem for a linearmodel for heterogeneousd-dimensional array data (d=1,2,3) organized inG known groups,and with identical tensor structured design matrixX across all groups. Specificallyn = \prod_i^d n_i is thenumber of observations in each group,Y_g then_1\times \cdots \times n_d response arrayfor groupg \in \{1,\ldots,G\}, andX an\times p design matrix, with tensor structure

X = \bigotimes_{i=1}^d X_i.

Ford =1,2,3,X_1,\ldots, X_d are the marginaln_i\times p_i design matrices (Kronecker components).Using the GLAM framework the model equation for groupg\in \{1,\ldots,G\} is expressed as

Y_g = \rho(X_d,\rho(X_{d-1},\ldots,\rho(X_1,B_g))) + E_g,

where\rho is the so called rotatedH-transfrom (seeCurrie et al., 2006),B_g for eachg is a (random)p_1\times\cdots\times p_d parameter arrayandE_g is an_1\times \cdots \times n_d error array.

This package solves the penalized soft maximin problem fromLund et al., 2020, given by

\min_{\beta}\frac{1}{\zeta}\log\bigg(\sum_{g=1}^G \exp(-\zeta \hat V_g(\beta))\bigg) + \lambda \Vert\beta\Vert_1, \quad \zeta > 0,\lambda \geq 0

for the setup described above. Note that

\hat V_g(\beta):=\frac{1}{n}(2\beta^\top X^\top vec(Y_g)-\beta^\top X^\top X\beta),

is the empirical explained variance fromMeinshausen and Buhlmann, 2015. SeeLund et al., 2020 for more details and references.

Ford=1,2,3, using only the marginal matricesX_1,X_2,\ldots (ford=1 there is only one marginal), the functionsoftmaximinsolves the soft maximin problem for a sequence of penalty parameters\lambda_{max}>\ldots >\lambda_{min}>0.

Two optimization algorithms are implemented, a non-monotoneproximal gradient (NPG) algorithm and a fast iterative soft thresholding algorithm (FISTA).We note that this package also solves the problem above with the penalty given by the SCADpenalty, using the multiple step adaptive lasso procedure to loop over the proximal algorithm.

Value

An object with S3 Class "SMMA".

Author(s)

Adam Lund

Maintainer: Adam Lund,adam.lund@math.ku.dk

References

Lund, A., S. W. Mogensen and N. R. Hansen (2020). Soft Maximin Estimation for Heterogeneous Array Data.Preprint.

Meinshausen, N and P. Buhlmann (2015). Maximin effects in inhomogeneous large-scale data.The Annals of Statistics. 43, 4, 1801-1830. url = https://doi.org/10.1214/15-AOS1325.

Currie, I. D., M. Durban, and P. H. C. Eilers (2006). Generalized lineararray models with applications to multidimensional smoothing.Journal of the Royal Statistical Society. Series B. 68, 259-280. url = http://dx.doi.org/10.1111/j.1467-9868.2006.00543.x.

Examples

##size of examplen1 <- 65; n2 <- 26; n3 <- 13; p1 <- 13; p2 <- 5; p3 <- 4##marginal design matrices (Kronecker components)X1 <- matrix(rnorm(n1 * p1), n1, p1)X2 <- matrix(rnorm(n2 * p2), n2, p2)X3 <- matrix(rnorm(n3 * p3), n3, p3)X <- list(X1, X2, X3)component <- rbinom(p1 * p2 * p3, 1, 0.1)Beta1 <- array(rnorm(p1 * p2 * p3, 0, 0.1) + component, c(p1 , p2, p3))mu1 <- RH(X3, RH(X2, RH(X1, Beta1)))Y1 <- array(rnorm(n1 * n2 * n3), dim = c(n1, n2, n3)) + mu1Beta2 <- array(rnorm(p1 * p2 * p3, 0, 0.1) + component, c(p1 , p2, p3))mu2 <- RH(X3, RH(X2, RH(X1, Beta2)))Y2 <- array(rnorm(n1 * n2 * n3), dim = c(n1, n2, n3)) + mu2Beta3 <- array(rnorm(p1 * p2 * p3, 0, 0.1) + component, c(p1 , p2, p3))mu3 <- RH(X3, RH(X2, RH(X1, Beta3)))Y3 <- array(rnorm(n1 * n2 * n3), dim = c(n1, n2, n3)) + mu3Beta4 <- array(rnorm(p1 * p2 * p3, 0, 0.1) + component, c(p1 , p2, p3))mu4 <- RH(X3, RH(X2, RH(X1, Beta4)))Y4 <- array(rnorm(n1 * n2 * n3), dim = c(n1, n2, n3)) + mu4Beta5 <- array(rnorm(p1 * p2 * p3, 0, 0.1) + component, c(p1 , p2, p3))mu5 <- RH(X3, RH(X2, RH(X1, Beta5)))Y5 <- array(rnorm(n1 * n2 * n3), dim = c(n1, n2, n3)) + mu5Y <- array(NA, c(dim(Y1), 5))Y[,,, 1] <- Y1; Y[,,, 2] <- Y2; Y[,,, 3] <- Y3; Y[,,, 4] <- Y4; Y[,,, 5] <- Y5;fit <- softmaximin(X, Y, zeta = 10, penalty = "lasso", alg = "npg")Betafit <- fit$coefmodelno <- 15m <- min(Betafit[ , modelno], c(component))M <- max(Betafit[ , modelno], c(component))plot(c(component), type="l", ylim = c(m, M))lines(Betafit[ , modelno], col = "red")

Make Prediction From a SMMA Object

Description

Given new covariate data this function computes the linear predictorsbased on the estimated model coefficients in an object produced by the functionsoftmaximin. Note that thedata can be supplied in two different formats: i) as an' \times p matrix (p is the number of modelcoefficients andn' is the number of new data points) or ii) as a list of two or three matrices each ofsizen_i' \times p_i, i = 1, 2, 3 (n_i' is the number of new marginal data points in theith dimension).

Usage

## S3 method for class 'SMMA'predict(object, x = NULL, X = NULL, ...)

Arguments

Value

A list of lengthnlambda containing the linear predictors for each model. Ifnew covariate data is supplied in onen' \times p matrixx eachitem is a vector of lengthn'. If the data is supplied as a list ofmatrices each of sizen'_{i} \times p_i, each item is an array of sizen'_1 \times \cdots \times n'_d, withd\in \{1,2,3\}.

Author(s)

Adam Lund

Examples

##size of examplen1 <- 65; n2 <- 26; n3 <- 13; p1 <- 13; p2 <- 5; p3 <- 4##marginal design matrices (Kronecker components)X1 <- matrix(rnorm(n1 * p1, 0, 0.5), n1, p1)X2 <- matrix(rnorm(n2 * p2, 0, 0.5), n2, p2)X3 <- matrix(rnorm(n3 * p3, 0, 0.5), n3, p3)X <- list(X1, X2, X3)component <- rbinom(p1 * p2 * p3, 1, 0.1)Beta1 <- array(rnorm(p1 * p2 * p3, 0, .1) + component, c(p1 , p2, p3))Beta2 <- array(rnorm(p1 * p2 * p3, 0, .1) + component, c(p1 , p2, p3))mu1 <- RH(X3, RH(X2, RH(X1, Beta1)))mu2 <- RH(X3, RH(X2, RH(X1, Beta2)))Y1 <- array(rnorm(n1 * n2 * n3, mu1), dim = c(n1, n2, n3))Y2 <- array(rnorm(n1 * n2 * n3, mu2), dim = c(n1, n2, n3))Y <- array(NA, c(dim(Y1), 2))Y[,,, 1] <- Y1; Y[,,, 2] <- Y2;fit <- softmaximin(X, Y, zeta = 10, penalty = "lasso", alg = "npg")##new data in matrix formx <- matrix(rnorm(p1 * p2 * p3), nrow = 1)predict(fit, x = x)[[15]]##new data in tensor component formX1 <- matrix(rnorm(p1), nrow = 1)X2 <- matrix(rnorm(p2), nrow = 1)X3 <- matrix(rnorm(p3), nrow = 1)predict(fit, X = list(X1, X2, X3))[[15]]

Print Function for objects of Class SMMA

Description

This function will print some information about the SMMA object.

Usage

## S3 method for class 'SMMA'print(x, ...)

Arguments

Author(s)

Adam Lund

Examples

 ##size of example n1 <- 65; n2 <- 26; n3 <- 13; p1 <- 13; p2 <- 5; p3 <- 4##marginal design matrices (Kronecker components)X1 <- matrix(rnorm(n1 * p1, 0, 0.5), n1, p1) X2 <- matrix(rnorm(n2 * p2, 0, 0.5), n2, p2) X3 <- matrix(rnorm(n3 * p3, 0, 0.5), n3, p3) X <- list(X1, X2, X3)component <- rbinom(p1 * p2 * p3, 1, 0.1) Beta1 <- array(rnorm(p1 * p2 * p3, 0, .1) + component, c(p1 , p2, p3))Beta2 <- array(rnorm(p1 * p2 * p3, 0, .1) + component, c(p1 , p2, p3))mu1 <- RH(X3, RH(X2, RH(X1, Beta1)))mu2 <- RH(X3, RH(X2, RH(X1, Beta2)))Y1 <- array(rnorm(n1 * n2 * n3, mu1), dim = c(n1, n2, n3))Y2 <- array(rnorm(n1 * n2 * n3, mu2), dim = c(n1, n2, n3))Y <- array(NA, c(dim(Y1), 2))Y[,,, 1] <- Y1; Y[,,, 2] <- Y2;fit <- softmaximin(X, Y, zeta = 10, penalty = "lasso", alg = "npg")fit