mNIX_marg: Marginal posterior of the mNIX hyperparameters as a...
In mlysy/losmix: Inference for Gaussian Location-Scale Mixed-Effects Models
Description
Usage
Arguments
Details
Value
Examples
Description
Marginal posterior of the mNIX hyperparameters as a TMB::MakeADFun object.
Usage
1
mnix_marg(id, y, X)
Arguments
id
Subject identifiers. A vector of length n consisting of nsub <= n distinct elements. If missing default to rep(1, n), i.e., only one subject.
y
Response vector for all subjects. A vector of length n.
X
Covariate matrix from all subjects. A matrix of size n x p.
Details
In particular, the return list contains elements fn, gr, and he corresponding to the negative log-posterior and its first and second derivatives, computed analytically in C++ via automatic differentiation (AD).
To evaluate these functions, the mNIX hyperparameters Phi = (lambda, Omega, nu, tau) must first be converted to an unconstrained vector, namely
1
where link{log_chol} is the log-Cholesky decomposition.
The default prior is
1
which is equivalent to
1
mnix_marg can also be used to simulate random effects beta_i, sigma_i) for each subject, using the simulate() method of the TMB::MakeADFun return object (see Examples).
Value
A list as returned by TMB::MakeADFun corresponding to the negative log-posterior of the marginal distribution of the mNIX hyperparameters (see Details).
Examples
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
# simulate data
set.seed(1234) # reproducible results
p <- 3 # number of covariates
nSub <- 50 # number of subjects
N <- sample(20:50, nSub, replace = TRUE) # observations per subject
# hyperparameters
lambda <- rnorm(p)
Omega <- diag(p)
nu <- 5 * runif(1) + 3
tau <- rexp(1)
# parameters
Theta <- mnix_sim(nSub, lambda = lambda, Omega = Omega, nu = nu, tau = tau)
# data
# design matrices
X <- lapply(N, function(n) matrix(rnorm(n*p), n, p))
# response vectors
y <- lapply(1:nSub, function(ii) {
rnorm(N[ii], mean = c(X[[ii]] %*% Theta$beta[ii,]), sd = Theta$sigma[ii])
})
# convert lists to single matrix/vector
X <- do.call(rbind, X)
y <- do.call(c, y)
# subject identifiers
id <- rep(1:nSub, times = N)
# marginal posterior distribution
nlp <- mnix_marg(id = id, y = y, X = X)
# mode-finding using quasi-Newton algorithm in nlm
# objective function
# note that gradient information is included via the "attribute",
# to increase computational efficiency when function and gradient
# share calculations.
# this is the case for TMB functions, as calling nlp$gr() without
# arguments reuses calculations from the last call to nlp$fn.
ofun <- function(par) {
out <- nlp$fn(par)
# include gradient information via 'attribute'
attr(out, "gradient") <- nlp$gr()
out
}
opt <- nlm(p = nlp$par, # starting value
f = ofun) # objective function
# gradient at the mode relative to its absolute value
nlp$gr(opt$estimate)/abs(opt$estimate)
## Approximate Bayesian Inference ##
#
# On transformed scale:
#
# - posterior mean is mode of nlp
# - posterior variance is inverse of hessian of nlp
#
# On original scale:
#
# - simulate from multivariate normal on the transformed scale
# - convert output to original scale
# transformed scale simulation
npost <- 1e4 # number of posterior draws
psi_mean <- opt$estimate # posterior mean
psi_var <- solveV(nlp$he(opt$estimate))
Psi_post <- rmvn(n = npost, mu = psi_mean, Sigma = psi_var)
# convert to original scale
Phi_post <- ivec_phi(Psi_post)
# histograms of posterior distributions
# true hyperparameter values are vertical lines in red
# format data for plotting
iOmega <- cbind(rbind(1:p, 1:p), combn(1:p,2)) # unique elements of Omega
Phi_plot <- cbind(Phi_post$lambda, # lambda
apply(iOmega, 2,
function(ii) Phi_post$Omega[ii[1],ii[2],]), # Omega
Phi_post$nu, # nu
Phi_post$tau) # tau
# hyperparameter names
phi_names <- c(paste0("lambda[", 1:p, "]"),
paste0("Omega[", iOmega[1,], iOmega[2,], "]"),
"nu", "tau")
# true values
phi_true <- c(lambda, Omega[t(iOmega)], nu, tau)
# create plot
par(mfrow = c(3,4), mar = c(2,2,4,.5))
for(ii in 1:ncol(Phi_plot)) {
# approximate posterior
hist(Phi_plot[,ii], breaks = 40, xlab = "", ylab = "",
main = parse(text = paste0("hat(p)(", phi_names[ii],
"*\" | \"*bold(Y),bold(X))")))
# true parameter value
abline(v = phi_true[ii], col = "red", lwd = 2)
}
# legend
plot(0, type = "n", xlim = c(0,1), ylim = c(0,1),
xlab = "", ylab = "", axes = FALSE)
legend("bottom", inset = .05,
legend = c("Posterior Distribution", "True Hyperparameter Value"),
lwd = c(NA, 2), pch = c(22, NA), seg.len = 1.5,
col = c("black", "red"), bg = c("white", NA), cex = .85)
## Inference for Random Effects ##
iSub <- sample(nSub, 1) # pick a subject at random
# data for subject i
Xi <- X[id == iSub,]
yi <- y[id == iSub]
# true parameter values for subject i
thetai_true <- c(Theta$beta[iSub,], sigma = Theta$sigma[iSub])
# method 1: using mnix_sim
set.seed(1)
Thetai_post <- mnix_sim(n = npost, y = yi, X = Xi,
lambda = Phi_post$lambda, Omega = Phi_post$Omega,
nu = Phi_post$nu, tau = Phi_post$tau)
# method 2: using mnix_marg
nlpi <- mnix_marg(y = yi, X = Xi) # posterior distribution given just subject i
set.seed(1)
Thetai_post2 <- apply(Psi_post, 1, function(psi) nlpi$simulate(psi))
# reshape data
Thetai_post2 <- unlist_bind(Thetai_post2,
name = c("beta", "sigma"),
bind = c(cbind, c))
Thetai_post2$beta <- t(Thetai_post2$beta) # beta needs to be transposed
# check that they are the same
all.equal(Thetai_post, Thetai_post2)
# format data for plotting
Thetai_plot <- cbind(Thetai_post$beta, # beta
Thetai_post$sigma) # sigma
# parameter names
thetai_names <- c(paste0("beta[i", 1:p, "]"), "sigma[i]")
# true parameter values for subject i
thetai_true <- c(Theta$beta[iSub,], sigma = Theta$sigma[iSub])
# create plot
par(mfrow = c(2,2), mar = c(2,2,4,.5))
for(ii in 1:ncol(Thetai_plot)) {
# approximate posterior
hist(Thetai_plot[,ii], breaks = 40, xlab = "", ylab = "",
main = parse(text = paste0("hat(p)(", thetai_names[ii],
"*\" | \"*bold(Y),bold(X))")))
# true parameter value
abline(v = thetai_true[ii], col = "red", lwd = 2)
}
mlysy/losmix documentation built on Jan. 18, 2021, 5:56 a.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.
Description Usage Arguments Details Value Examples
Description
Marginal posterior of the mNIX hyperparameters as a TMB::MakeADFun object.
Usage
1 | mnix_marg(id, y, X)
|
Arguments
id |
Subject identifiers. A vector of length |
y |
Response vector for all subjects. A vector of length |
X |
Covariate matrix from all subjects. A matrix of size |
Details
In particular, the return list contains elements fn, gr, and he corresponding to the negative log-posterior and its first and second derivatives, computed analytically in C++ via automatic differentiation (AD).
To evaluate these functions, the mNIX hyperparameters Phi = (lambda, Omega, nu, tau) must first be converted to an unconstrained vector, namely
1 |
where link{log_chol} is the log-Cholesky decomposition.
The default prior is
1 |
which is equivalent to
1 |
mnix_marg can also be used to simulate random effects beta_i, sigma_i) for each subject, using the simulate() method of the TMB::MakeADFun return object (see Examples).
Value
A list as returned by TMB::MakeADFun corresponding to the negative log-posterior of the marginal distribution of the mNIX hyperparameters (see Details).
Examples
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 | # simulate data
set.seed(1234) # reproducible results
p <- 3 # number of covariates
nSub <- 50 # number of subjects
N <- sample(20:50, nSub, replace = TRUE) # observations per subject
# hyperparameters
lambda <- rnorm(p)
Omega <- diag(p)
nu <- 5 * runif(1) + 3
tau <- rexp(1)
# parameters
Theta <- mnix_sim(nSub, lambda = lambda, Omega = Omega, nu = nu, tau = tau)
# data
# design matrices
X <- lapply(N, function(n) matrix(rnorm(n*p), n, p))
# response vectors
y <- lapply(1:nSub, function(ii) {
rnorm(N[ii], mean = c(X[[ii]] %*% Theta$beta[ii,]), sd = Theta$sigma[ii])
})
# convert lists to single matrix/vector
X <- do.call(rbind, X)
y <- do.call(c, y)
# subject identifiers
id <- rep(1:nSub, times = N)
# marginal posterior distribution
nlp <- mnix_marg(id = id, y = y, X = X)
# mode-finding using quasi-Newton algorithm in nlm
# objective function
# note that gradient information is included via the "attribute",
# to increase computational efficiency when function and gradient
# share calculations.
# this is the case for TMB functions, as calling nlp$gr() without
# arguments reuses calculations from the last call to nlp$fn.
ofun <- function(par) {
out <- nlp$fn(par)
# include gradient information via 'attribute'
attr(out, "gradient") <- nlp$gr()
out
}
opt <- nlm(p = nlp$par, # starting value
f = ofun) # objective function
# gradient at the mode relative to its absolute value
nlp$gr(opt$estimate)/abs(opt$estimate)
## Approximate Bayesian Inference ##
#
# On transformed scale:
#
# - posterior mean is mode of nlp
# - posterior variance is inverse of hessian of nlp
#
# On original scale:
#
# - simulate from multivariate normal on the transformed scale
# - convert output to original scale
# transformed scale simulation
npost <- 1e4 # number of posterior draws
psi_mean <- opt$estimate # posterior mean
psi_var <- solveV(nlp$he(opt$estimate))
Psi_post <- rmvn(n = npost, mu = psi_mean, Sigma = psi_var)
# convert to original scale
Phi_post <- ivec_phi(Psi_post)
# histograms of posterior distributions
# true hyperparameter values are vertical lines in red
# format data for plotting
iOmega <- cbind(rbind(1:p, 1:p), combn(1:p,2)) # unique elements of Omega
Phi_plot <- cbind(Phi_post$lambda, # lambda
apply(iOmega, 2,
function(ii) Phi_post$Omega[ii[1],ii[2],]), # Omega
Phi_post$nu, # nu
Phi_post$tau) # tau
# hyperparameter names
phi_names <- c(paste0("lambda[", 1:p, "]"),
paste0("Omega[", iOmega[1,], iOmega[2,], "]"),
"nu", "tau")
# true values
phi_true <- c(lambda, Omega[t(iOmega)], nu, tau)
# create plot
par(mfrow = c(3,4), mar = c(2,2,4,.5))
for(ii in 1:ncol(Phi_plot)) {
# approximate posterior
hist(Phi_plot[,ii], breaks = 40, xlab = "", ylab = "",
main = parse(text = paste0("hat(p)(", phi_names[ii],
"*\" | \"*bold(Y),bold(X))")))
# true parameter value
abline(v = phi_true[ii], col = "red", lwd = 2)
}
# legend
plot(0, type = "n", xlim = c(0,1), ylim = c(0,1),
xlab = "", ylab = "", axes = FALSE)
legend("bottom", inset = .05,
legend = c("Posterior Distribution", "True Hyperparameter Value"),
lwd = c(NA, 2), pch = c(22, NA), seg.len = 1.5,
col = c("black", "red"), bg = c("white", NA), cex = .85)
## Inference for Random Effects ##
iSub <- sample(nSub, 1) # pick a subject at random
# data for subject i
Xi <- X[id == iSub,]
yi <- y[id == iSub]
# true parameter values for subject i
thetai_true <- c(Theta$beta[iSub,], sigma = Theta$sigma[iSub])
# method 1: using mnix_sim
set.seed(1)
Thetai_post <- mnix_sim(n = npost, y = yi, X = Xi,
lambda = Phi_post$lambda, Omega = Phi_post$Omega,
nu = Phi_post$nu, tau = Phi_post$tau)
# method 2: using mnix_marg
nlpi <- mnix_marg(y = yi, X = Xi) # posterior distribution given just subject i
set.seed(1)
Thetai_post2 <- apply(Psi_post, 1, function(psi) nlpi$simulate(psi))
# reshape data
Thetai_post2 <- unlist_bind(Thetai_post2,
name = c("beta", "sigma"),
bind = c(cbind, c))
Thetai_post2$beta <- t(Thetai_post2$beta) # beta needs to be transposed
# check that they are the same
all.equal(Thetai_post, Thetai_post2)
# format data for plotting
Thetai_plot <- cbind(Thetai_post$beta, # beta
Thetai_post$sigma) # sigma
# parameter names
thetai_names <- c(paste0("beta[i", 1:p, "]"), "sigma[i]")
# true parameter values for subject i
thetai_true <- c(Theta$beta[iSub,], sigma = Theta$sigma[iSub])
# create plot
par(mfrow = c(2,2), mar = c(2,2,4,.5))
for(ii in 1:ncol(Thetai_plot)) {
# approximate posterior
hist(Thetai_plot[,ii], breaks = 40, xlab = "", ylab = "",
main = parse(text = paste0("hat(p)(", thetai_names[ii],
"*\" | \"*bold(Y),bold(X))")))
# true parameter value
abline(v = thetai_true[ii], col = "red", lwd = 2)
}
|
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.