Skip to content

Commit ac87ec4

Browse files
author
Steve Walker
committed
JSS version
1 parent 23fec45 commit ac87ec4

File tree

4 files changed

+167
-1
lines changed

4 files changed

+167
-1
lines changed

‎NAMESPACE

+1
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ export(mkZt)
1414
export(mkZtSection)
1515
export(pirls)
1616
export(pls)
17+
export(plsJSS)
1718
export(plsform)
1819
export(rLmer)
1920
importFrom(Matrix,Cholesky)

‎R/JSS.R

+121
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,121 @@
1+
##' Linear mixed model deviance function as it
2+
##' appears in the pseudocode of the JSS article
3+
##'
4+
##' A pure \code{R} implementation of the
5+
##' penalized least squares (PLS) approach for computing
6+
##' linear mixed model deviances. The purpose
7+
##' is to clarify how PLS works without having
8+
##' to read through C++ code, and as a sandbox for
9+
##' trying out modifications to PLS.
10+
##'
11+
##' @param X fixed effects model matrix
12+
##' @param y response
13+
##' @param Zt transpose of the sparse model matrix for the random effects
14+
##' @param Lambdat upper triangular sparse Cholesky factor of the
15+
##' relative covariance matrix of the random effects
16+
##' @param mapping a function that takes a value of \code{theta} and produces
17+
##' the non-zero elements of \code{Lambdat}. The structure of \code{Lambdat}
18+
##' cannot change, only the numerical values
19+
##' @param weights prior weights
20+
##' @param offset offset
21+
##' @param REML calculate REML deviance?
22+
##' @param ... additional arguments
23+
##' @keywords models
24+
##'
25+
##' @return a function that evaluates the deviance or REML criterion
26+
##' @export
27+
plsJSS <- function(X, y, Zt, Lambdat, mapping, weights,
28+
offset = numeric(n), REML = TRUE, ...)
29+
{
30+
# SW: how to test for sparse matrices, without specifying the specific class?
31+
stopifnot(is.matrix(X)) # is.matrix(Zt), is.matrix(Lambdat))
32+
n <- length(y); p <- ncol(X); q <- nrow(Zt)
33+
stopifnot(nrow(X) == n, ncol(Zt) == n,
34+
nrow(Lambdat) == q, ncol(Lambdat) == q)
35+
# calculate weighted products
36+
sqrtW <- if (missing(weights)) Diagonal(n=n) else Diagonal(x=sqrt(as.numeric(weights)))
37+
WX <- sqrtW %*% X
38+
Wy <- sqrtW %*% y
39+
ZtW <- Zt %*% sqrtW
40+
XtWX <- crossprod(WX)
41+
XtWy <- crossprod(WX, Wy)
42+
ZtWX <- ZtW %*% WX
43+
ZtWy <- ZtW %*% Wy
44+
rm(WX,Wy)
45+
local({ # mutable values stored in local environment
46+
b <- numeric(q) # conditional mode of random effects
47+
beta <- numeric(p) # conditional estimate of fixed-effects
48+
cu <- numeric(q) # intermediate solution
49+
RXtRX <- XtWX # down-dated XtWX
50+
L <- Cholesky(tcrossprod(Lambdat %*% ZtW), LDL = FALSE, Imult=1)
51+
Lambdat <- Lambdat # stored here b/c x slot will be updated
52+
mu <- numeric(n) # conditional mean of response
53+
RZX <- matrix(0,nrow=q,ncol=p) # intermediate matrix in solution
54+
u <- numeric(q) # conditional mode of spherical random effects
55+
degFree <- as.numeric(n) # degrees of freedom (depends on REML)
56+
if(REML) degFree <- degFree - as.numeric(p)
57+
function(theta) {
58+
59+
##################################################
60+
# Step I: update covariance parameters
61+
##################################################
62+
# update relative covariance factor
63+
# by placing the new values of theta
64+
# in the appropriate positions
65+
Lambdat@x[] <<- mapping(theta)
66+
# update random-effects Cholesky factor
67+
L <<- update(L, Lambdat %*% ZtW, mult = 1)
68+
69+
70+
##################################################
71+
# Step II: solve normal equations
72+
##################################################
73+
# solve eqn. ??
74+
cu[] <<- as.vector(solve(L, solve(L, Lambdat %*% ZtWy,
75+
system="P"), system="L"))
76+
# solve eqn. ??
77+
RZX[] <<- as.vector(solve(L, solve(L, Lambdat %*% ZtWX,
78+
system="P"), system="L"))
79+
# downdate XtWX and form Cholesky
80+
# factor (eqn. ??)
81+
RXtRX <<- as(XtWX - crossprod(RZX), "dpoMatrix")
82+
# conditional estimate of fixed-effects
83+
# coefficients (solve eqn. ??)
84+
beta[] <<- as.vector(solve(RXtRX, XtWy - crossprod(RZX, cu)))
85+
# conditional mode of the spherical
86+
# random-effects coefficients (eqn. ??)
87+
u[] <<- as.vector(solve(L, solve(L, cu - RZX %*% beta,
88+
system = "Lt"), system="Pt"))
89+
# update conditional model of the
90+
# non-spherical random-effects
91+
# coefficients
92+
b[] <<- as.vector(crossprod(Lambdat,u))
93+
94+
95+
##################################################
96+
# Step III: update linear predictor and residuals
97+
##################################################
98+
# update linear predictor
99+
mu[] <<- as.vector(crossprod(Zt,b) + X %*% beta + offset)
100+
# weighted residuals
101+
wtres <- sqrtW*(y-mu)
102+
# penalized, weighted residual
103+
# sum-of-squares
104+
105+
106+
##################################################
107+
# Step IV: compute profiled deviance
108+
##################################################
109+
pwrss <- sum(wtres^2) + sum(u^2)
110+
# log determinant (depends on
111+
# whether REML or ML is used)
112+
logDet <- 2*determinant(L, logarithm = TRUE)$modulus
113+
if (REML) logDet <- logDet + determinant(RXtRX,
114+
logarithm = TRUE)$modulus
115+
attributes(logDet) <- NULL
116+
# profiled deviance or REML criterion
117+
profDev <- logDet + degFree*(1 + log(2*pi*pwrss) - log(degFree))
118+
return(profDev)
119+
}
120+
})
121+
}

‎man/plsJSS.Rd

+44
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
\name{plsJSS}
2+
\alias{plsJSS}
3+
\title{Linear mixed model deviance function as it
4+
appears in the pseudocode of the JSS article}
5+
\usage{
6+
plsJSS(X, y, Zt, Lambdat, mapping, weights, offset = numeric(n),
7+
REML = TRUE, ...)
8+
}
9+
\arguments{
10+
\item{X}{fixed effects model matrix}
11+
12+
\item{y}{response}
13+
14+
\item{Zt}{transpose of the sparse model matrix for the
15+
random effects}
16+
17+
\item{Lambdat}{upper triangular sparse Cholesky factor of
18+
the relative covariance matrix of the random effects}
19+
20+
\item{mapping}{a function that takes a value of
21+
\code{theta} and produces the non-zero elements of
22+
\code{Lambdat}. The structure of \code{Lambdat} cannot
23+
change, only the numerical values}
24+
25+
\item{weights}{prior weights}
26+
27+
\item{offset}{offset}
28+
29+
\item{REML}{calculate REML deviance?}
30+
31+
\item{...}{additional arguments}
32+
}
33+
\value{
34+
a function that evaluates the deviance or REML criterion
35+
}
36+
\description{
37+
A pure \code{R} implementation of the penalized least
38+
squares (PLS) approach for computing linear mixed model
39+
deviances. The purpose is to clarify how PLS works without
40+
having to read through C++ code, and as a sandbox for
41+
trying out modifications to PLS.
42+
}
43+
\keyword{models}
44+

‎tests/cbpp.R

+1-1
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,7 @@ data(cbpp, package = 'lme4')
1010
ll <- plsform(form, data = cbpp, family = binomial)
1111
devf <- do.call(pirls, c(ll, list(family=binomial)))
1212
devf(c(1,0,0,0,0))
13-
bobyqa(c(1,0,0,0,0), devf)
13+
(opt <- bobyqa(c(1,0,0,0,0), devf))
1414

1515
# get initial values
1616
gm1 <- glm(nobars(form), binomial, cbpp)

0 commit comments

Comments
 (0)