General model
Consider an abundance table \(Y\in
\mathbb{N}^{n\times p}\), with missing values. Assume that three
other tables are available: \(R\in
\mathbb{R}^{n\times K_1}\), containing covariates about the rows
of \(Y\) (geographical information for
instance), \(C\in \mathbb{R}^{p\times
K_2}\), containing covariates about the columns of \(Y\) (yearly meteorological indices), and
\(E\in \mathbb{R}^{np\times K_3}\)
containing co- variates about the row-column pairs (yearly
meteorological information at the sites’ scale, yearly economical
indices of the sites’ country, etc.). We denote by \(R_i\) (resp. \(C_i\), \(E_i\)) the \(i\)-th row of \(R\) (resp of \(C\), \(E\)).
The model implemented in the lori package is the following log-linear
model: \[Y_{i,j}\sim\mathcal{P}(\exp(X_{i,j})), \quad
X_{i,j} = \alpha_i + \beta_j + \sum_{k=1}^{K_1+K_2}U^{i,j}_k\epsilon_k +
\Theta_{i,j},
\] where \(U \in
\mathbb{R}^{np\times(K1+K2+K3)}\) denotes the covariate matrix
concatenating all the available covariate information, and whose \((j-1)n+i\)-th row, for \(1\leq i\leq n\) \(1\leq j\leq p\) are defined by: \(U_{(j-1)n+i} = (R_i, C_j, E_{j-1}n+i)\).
Note that the function covmat is implemented to build \(U\) from any combination of \(R\), \(C\)
and \(E\).
Simulated data
We start by generating simulated data as follows
## covariates
n <- 10 # number of rows
p <- 6 # number of columns
K1 <- 2 # number of row covariates
K2 <- 2 # number of column covariates
K3 <- 3 # number of (rowxcolumn) covariates
R <- matrix(rnorm(n*K1), nrow=n) # matrix of row covariates
C <- matrix(rnorm(p*K2), nrow=p) # matrix of column covariates
E <- matrix(rnorm(n*p*K3), nrow=n*p) # matrix of (rowxcolumn) covariates
U <- covmat(n, p, R, C, E)
U <- scale(U)
Note that the function also supports creating U from any singleton or
pair of the matrices \((R,C,E)\),
provided that the arguments are in the correct order. In other words,
all the commands below may also be used:
U <- covmat(n, p, R)
U <- covmat(n, p, C)
U <- covmat(n, p, R, C)
U <- covmat(n, p, R=R, E=E)
U <- covmat(n, p, C=C, E=E)
Then we simulate the parameters of the models.
## parameters
alpha0 <- rep(0, n)
alpha0[1:6] <- 1
beta0 <- rep(0, p)
beta0[1:4] <- 1
epsilon0 <- rep(0, K1+K2+K3)
epsilon0[5:6] <- 0.2
r <- 2 #rank of interaction matrix theta0
theta0 <- 0.1*matrix(rnorm(n*r), nrow=n)%*%diag(r:1)%*%matrix(rnorm(p*r), nrow=r)
theta0 <- sweep(theta0, 2, colMeans(theta0))
theta0 <- sweep(theta0, 1, rowMeans(theta0))
Finally we build the matrix of log-intensities, sample corresponding
Poisson data, and add MCAR missing values.
## construct x0
x0 <- matrix(rep(alpha0, p), nrow=n) #row effects
x0 <- x0 + matrix(rep(beta0, each=n), nrow=n) #add column effects
x0 <- x0 + matrix(U%*% epsilon0, nrow=n) #add cov effects
x0 <- x0 + theta0 #add interactions
## sample count data y
y0 <- matrix(rpois(n*p, lambda = c(exp(x0))), nrow = n)
## add missing values
p_miss <- 0.2
y <- y0
y[sample(1:(n*p), round(p_miss*n*p))] <- NA
# Estimation and imputation
The purpose of the lori package is dual: to estimate the parameters
\(\alpha\), \(\beta\), \(\epsilon\), \(\Theta\), and to impute the count table
\(Y\). The corresponding input are the
incomplete count data \(Y\) and
(optionally) the covariate matrix \(U\). The estimation procedure is: \[(\hat\alpha, \hat\beta, \hat\epsilon, \hat
\Theta)\in \operatorname{argmin} \mathcal{L}(Y; \alpha, \beta, \epsilon,
\Theta) + \lambda_1\left\|\Theta \right\|_* +
\lambda_2\left(\left\|\alpha \right\|_1+\left\|\beta
\right\|_1+\left\|\epsilon \right\|_1\right),
\] where \(\mathcal{L}(Y; \alpha,
\beta, \epsilon, \Theta)\) is the Poisson negative
log-likelihood. With \(\Omega\) the
mask encoding missing values (\(\Omega_{ij}
=1\) if entry \((i,j)\) is
observed and \(0\) otherwise): \[
\sum_{(i,j)}\Omega_{i,j}[-Y_{i,j}(\alpha_i+\beta_j +
\sum_{k=1}^{K}\epsilon_{k}U(i,j)_k + \Theta_{i,j})
+ \exp(\alpha_i+\beta_j + \sum_{k=1}^{K}\epsilon_{k}U(i,j)_k +
\Theta_{i,j})],
\] The intuition is that the negative log-likelihood is minimized
(the parameters should fit the data as much as possible) with additional
regularization terms which constrain the parameters and perform model
selection automatically. The first penalty term, \(\norm{\Theta}[*]\) corresponds to the
nuclear norm of \(\Theta\) and induces
low-rank solutions: this is interpreted as assuming that a few latent
factors summarize the interactions. The second term, \(\left\|\alpha \right\|_1+\left\|\beta
\right\|_1+\left\|\epsilon \right\|_1\), corresponds to the sum
of the \(\ell_1\) norms of the vectors
\(\alpha\), \(\beta\) and \(\epsilon\): \[\left\|\alpha \right\|_1+\left\|\beta
\right\|_1+\left\|\epsilon \right\|_1 = \sum_{i=1}^{m_1}|\alpha_i| +
\sum_{j=1}^{m_2}|\beta_j| + \sum_{k=1}^{K}|\epsilon_k|.\] This
penalty induces sparse solutions (vectors \(\alpha\), \(\beta\) and \(\epsilon\) containing many zeros), meaning
that not all rows, columns and covariates have an effect on the counts.
The trade-off between the data-fitting term and the penalties is
controlled by the regularization parameters \(\lambda_1\) and \(\lambda_2\). As \(\lambda_1\) increases, the rank of the
solution \(\hat\Theta\) decreases (the
number of latent factors decreases). As \(\lambda_2\) increases, the number of
nonzero values in \(\hat\alpha\), \(\hat\beta\) and \(\hat\epsilon\) decreases (the number of
active rows, columns and covariates decreases).
In R, the estimation is done as follows, for some predefined
regularization parameters \(\lambda_1\)
and \(\lambda_2\):
## lori estimation
lambda1 <- 0.1
lambda2 <- 0.1
m <- sum(!is.na(y))
t <- Sys.time()
res.lori <- lori(y, U, 0.1, 0.1)
t <- Sys.time()-t
However, the quality of the estimation is highly dependent on the
choice of \(\lambda_1\) and \(\lambda_2\), and thus we provide a
cross-validation function to select them automatically:
res.cv <- cv.lori(y, U, trace.it = F, len=5)
res.lori <- lori(y, U, res.cv$lambda1, res.cv$lambda2)
The estimated parameters are accessed as follows
res.lori$alpha
#> [1] 0.13017010 0.42375452 0.00000000 0.32546414 0.38532507 -0.01740698
#> [7] -0.49502015 -0.08910152 -0.24525029 -0.40457170
res.lori$beta
#> [1] 0.1015438 0.1393790 0.3400212 0.4506407 -0.9235591 -0.2277350
res.lori$espilon
#> NULL
res.lori$theta
#> [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.56858957 -0.148792558 -0.029314118 0.074562971 -0.030270667
#> [2,] 0.78081972 -0.044675990 0.439917938 -0.328623857 -0.101531444
#> [3,] -0.75828423 0.255054213 0.209386276 -0.252298239 0.019104295
#> [4,] -1.01427258 0.005845634 -0.728405490 0.567768208 0.151488118
#> [5,] 0.25273257 -0.111536099 -0.149571762 0.155706536 0.003595749
#> [6,] -0.78660964 0.387798332 0.587791968 -0.594369899 -0.026459027
#> [7,] 0.09343521 -0.178513505 -0.468173273 0.428174794 0.052887602
#> [8,] -0.11092376 0.055522130 0.085404073 -0.086074034 -0.004045380
#> [9,] 0.46676145 -0.112400707 0.005243828 0.034901749 -0.028509395
#> [10,] 0.50775169 -0.108301451 0.047720560 0.000251772 -0.036259850
#> [,6]
#> [1,] -0.43477520
#> [2,] -0.74590636
#> [3,] 0.52703769
#> [4,] 1.01757611
#> [5,] -0.15092699
#> [6,] 0.43184826
#> [7,] 0.07218918
#> [8,] 0.06011697
#> [9,] -0.36599693
#> [10,] -0.41116272
The function also returns imputed values for the missing entries in
\(\Y\), based on these point estimates.
For each missing entry, the predicted value is given by \[\hat{Y}_{i,j} = \exp\left(\hat\alpha_i +
\hat\beta_j + \sum_{k=1}^{K_1+K_2}U_{(j-1)n+i,k}\hat\epsilon_k +
\hat\Theta_{i,j}\right).
\] The imputed data set \(\hatY\), such that \(\hat{Y}_{i,j} = Y_{i,j}\) if \(Y_{i,j}\) is observed, and \(\hat{Y}_{i,j}\) given by the equation above
otherwise, is accessed as follows:
res.lori$imputed
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 9 3 2 5 2 1
#> [2,] 11 10 14 2 1 1
#> [3,] 4 5 9 8 0 2
#> [4,] 1 4 4 9 4 9
#> [5,] 14 4 4 12 1 1
#> [6,] 4 5 9 4 2 5
#> [7,] 1 3 0 7 2 2
#> [8,] 0 7 2 3 0 0
#> [9,] 3 2 3 6 2 1
#> [10,] 4 2 2 2 0 1
