SCM

SCM Repository

[matrix] Annotation of /pkg/src/lmer.c
ViewVC logotype

Annotation of /pkg/src/lmer.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1162 - (view) (download) (as text)

1 : bates 415 #include "lmer.h"
2 : bates 808 #include <float.h>
3 : bates 255
4 : bates 1143 /* These should be moved to the R sources */
5 : maechler 534
6 :     /**
7 : bates 1143 * Symmetrize a matrix by copying the strict upper triangle into the
8 :     * lower triangle.
9 : maechler 534 *
10 : bates 1143 * @param a pointer to a matrix in Fortran storage mode
11 :     * @param nc number of columns (and rows and leading dimension) in the matrix
12 : maechler 534 *
13 : bates 1143 * @return a, symmetrized
14 : bates 275 */
15 : bates 1143 static double*
16 :     internal_symmetrize(double *a, int nc)
17 : bates 275 {
18 : bates 1143 int i, j;
19 : maechler 534
20 : bates 1143 for (i = 1; i < nc; i++)
21 :     for (j = 0; j < i; j++)
22 :     a[i + j*nc] = a[j + i*nc];
23 :     return a;
24 : bates 255 }
25 :    
26 : bates 742 /**
27 : bates 1143 * Simulate the Cholesky factor of a standardized Wishart variate with
28 :     * dimension p and df degrees of freedom.
29 : bates 742 *
30 : bates 1143 * @param df degrees of freedom
31 :     * @param p dimension of the Wishart distribution
32 :     * @param ans array of size p * p to hold the result
33 : bates 742 *
34 : bates 1143 * @return ans
35 : bates 742 */
36 : bates 1143 static double*
37 :     std_rWishart_factor(double df, int p, double ans[])
38 : bates 742 {
39 : bates 1143 int i, j, pp1 = p + 1;
40 : bates 742
41 : bates 1143 if (df < (double) p || p <= 0)
42 :     error("inconsistent degrees of freedom and dimension");
43 :     for (j = 0; j < p; j++) { /* jth column */
44 :     ans[j * pp1] = sqrt(rchisq(df - (double) j));
45 :     for (i = 0; i < j; i++) ans[i + j * p] = norm_rand();
46 : bates 742 }
47 : bates 1143 return ans;
48 : bates 742 }
49 :    
50 :    
51 : bates 1143 /* Internally used utilities */
52 :    
53 :     #define flag_not_factored(x) *LOGICAL(GET_SLOT(x, Matrix_statusSym)) = 0
54 :    
55 :     /**
56 :     * Create the diagonal blocks of the variance-covariance matrix of the
57 :     * random effects
58 :     *
59 :     * @param Linv - cholmod_sparse representation of L^{-1}
60 :     * @param nf - number of grouping factors
61 :     * @param Gp - group pointers
62 :     * @param nc - number of columns per factor
63 :     * @param bVar - list of 3-d arrays to be filled in
64 :     * @param uplo - "U" or "L" for upper or lower triangle
65 : bates 255 */
66 : bates 1143 static void
67 :     Linv_to_bVar(cholmod_sparse *Linv, const int Gp[], const int nc[],
68 :     SEXP bVar, const char uplo[])
69 : bates 255 {
70 : bates 1143 int *Lii = (int*)(Linv->i), *Lip = (int*)(Linv->p), i, nf = LENGTH(bVar);
71 :     double *Lix = (double*)(Linv->x), one[] = {1,0}, zero[] = {0,0};
72 :    
73 : bates 268 for (i = 0; i < nf; i++) {
74 : bates 1143 int *ind, j, nci = nc[i], maxnnz = 0;
75 :     int ncisqr = nci * nci, nlev = (Gp[i + 1] - Gp[i])/nci;
76 :     double *bVi = REAL(VECTOR_ELT(bVar, i)), *tmp;
77 : maechler 534
78 : bates 1143 AZERO(bVi, nlev * ncisqr);
79 :     for (j = 0; j < nlev; j++) {
80 :     int nzm = Lip[Gp[i] + (j + 1) * nci] - Lip[Gp[i] + j * nci];
81 :     if (nzm > maxnnz) maxnnz = nzm;
82 : bates 255 }
83 : bates 1143 ind = Calloc(maxnnz, int);
84 :     tmp = Calloc(maxnnz * nci, double);
85 :     for (j = 0; j < nlev; j++) {
86 :     int jj, k, kk;
87 :     int *ap = Lip + Gp[i] + j * nci;
88 :     int nr = ap[1] - ap[0];
89 : bates 255
90 : bates 1143 AZERO(tmp, maxnnz * nci);
91 :     Memcpy(ind, Lii + ap[0], nr);
92 :     Memcpy(tmp, Lix + ap[0], nr);
93 :     for (jj = 1; jj < nci; jj++) {
94 :     for (k = ap[jj]; k < ap[jj + 1]; k++) {
95 :     int aik = Lii[k];
96 :     for (kk = 0; kk < nr; kk++) {
97 :     if (aik == ind[kk]) {
98 :     tmp[kk + jj * maxnnz] = Lix[k];
99 :     aik = -1;
100 :     break;
101 :     }
102 :     }
103 :     if (aik >= 0) { /* did not find the row index */
104 :     ind[nr] = aik;
105 :     tmp[nr + jj * maxnnz] = Lix[k];
106 :     nr++;
107 :     }
108 :     }
109 : bates 901 }
110 : bates 1143 F77_CALL(dsyrk)(uplo, "T", &nci, &nr, one, tmp, &maxnnz,
111 :     zero, bVi + j * ncisqr, &nci);
112 : bates 901 }
113 : bates 1143 Free(ind); Free(tmp);
114 : bates 901 }
115 :     }
116 :    
117 : bates 1162 static void
118 :     internal_mer_bVar(SEXP x)
119 :     {
120 :     int q = LENGTH(GET_SLOT(x, Matrix_rZySym));
121 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
122 :     cholmod_sparse *eye = cholmod_speye(q, q, CHOLMOD_REAL, &c), *Linv;
123 :     int *Perm = (int *)(L->Perm), *iperm = Calloc(q, int), i;
124 :     /* create the inverse permutation */
125 :     for (i = 0; i < q; i++) iperm[Perm[i]] = i;
126 :     /* apply iperm to the identity matrix */
127 :     for (i = 0; i < q; i++) ((int*)(eye->i))[i] = iperm[i];
128 :     /* Create Linv */
129 :     Linv = cholmod_spsolve(CHOLMOD_L, L, eye, &c);
130 :     cholmod_free_sparse(&eye, &c);
131 :     Linv_to_bVar(Linv, INTEGER(GET_SLOT(x, Matrix_GpSym)),
132 :     INTEGER(GET_SLOT(x, Matrix_ncSym)),
133 :     GET_SLOT(x, Matrix_bVarSym), "U");
134 :     cholmod_free_sparse(&Linv, &c);
135 :     Free(L);
136 :     }
137 :    
138 : bates 1143 /**
139 :     * Evaluate the quadratic form in b defined by Omega
140 :     *
141 :     * @param b vector of random effects
142 :     * @param Omega - list of dpoMatrix objects defining the pattern for Omega
143 :     * @param nf - number of grouping factors
144 :     * @param Gp - group pointers
145 :     * @param nc - number of columns per factor
146 :     *
147 :     * @return
148 : bates 300 */
149 : bates 1143 static double
150 :     b_quadratic(const double b[], SEXP Omega, const int Gp[], const int nc[])
151 : bates 300 {
152 : bates 1143 int i, ione = 1, nf = LENGTH(Omega);
153 :     double ans = 0., one[] = {1.,0.};
154 : bates 300
155 : bates 415 for (i = 0; i < nf; i++) {
156 : bates 1143 int nci = nc[i], ntot = Gp[i + 1] - Gp[i];
157 :     int nlev = ntot/nci;
158 :     double *bcp = Memcpy(Calloc(ntot, double), b + Gp[i], ntot),
159 :     *omgf = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)), Matrix_xSym));
160 : bates 300
161 : bates 1143 F77_CALL(dtrmm)("L", "U", "N", "N", &nci, &nlev, one, omgf, &nci, bcp, &nci);
162 :     ans += F77_CALL(ddot)(&ntot, bcp, &ione, bcp, &ione);
163 :     Free(bcp);
164 : bates 274 }
165 : bates 1143 return ans;
166 : bates 274 }
167 :    
168 : maechler 534 /**
169 : bates 1143 * Calculate the length of the parameter vector (historically called "coef"
170 :     * even though these are not coefficients).
171 : maechler 534 *
172 : bates 1143 * @param nf number of factors
173 :     * @param nc number of columns in the model matrices for each factor
174 : maechler 534 *
175 : bates 1143 * @return total length of the coefficient vector
176 : bates 445 */
177 : bates 447 static R_INLINE
178 : bates 1143 int coef_length(int nf, const int nc[])
179 : bates 445 {
180 : bates 1143 int i, ans = 0;
181 :     for (i = 0; i < nf; i++) ans += (nc[i] * (nc[i] + 1))/2;
182 :     return ans;
183 : bates 445 }
184 : bates 291
185 : bates 1143 /**
186 :     * Find a variable of a given name in a given environment and check
187 :     * that its length and mode are correct.
188 :     *
189 :     * @param rho Environment in which to find the variable
190 :     * @param nm Name of the variable to find
191 :     * @param mode Desired mode
192 :     * @param len Desired length
193 :     *
194 :     * @return
195 : bates 255 */
196 : bates 1143 static
197 :     SEXP find_and_check(SEXP rho, SEXP nm, SEXPTYPE mode, int len)
198 :     {
199 :     SEXP ans;
200 :     if (R_NilValue == PROTECT(ans = findVarInFrame(rho, nm)))
201 :     error(_("environment `rho' must contain an object `%s'"),
202 :     CHAR(PRINTNAME(nm)));
203 :     if (TYPEOF(ans) != mode)
204 :     error(_("object `%s' of incorrect type"),
205 :     CHAR(PRINTNAME(nm)));
206 :     if (len && LENGTH(ans) != len)
207 :     error(_("object `%s' must be of length `%d'"),
208 :     CHAR(PRINTNAME(nm)), len);
209 :     UNPROTECT(1);
210 :     return ans;
211 : bates 255 }
212 :    
213 : bates 1143 /**
214 :     * Evaluate an expression in an environment, check that the length and
215 :     * mode are as expected and store the result.
216 :     *
217 :     * @param fcn expression to evaluate
218 :     * @param rho environment in which to evaluate it
219 :     * @param vv position to store the result
220 :     *
221 :     * @return vv with new contents
222 : bates 415 */
223 : bates 1143 static
224 :     SEXP eval_check_store(SEXP fcn, SEXP rho, SEXP vv)
225 : bates 415 {
226 : bates 1143 SEXP v = PROTECT(eval(fcn, rho));
227 :     if (TYPEOF(v) != TYPEOF(vv) || LENGTH(v) != LENGTH(vv))
228 :     error(_("fcn produced mode %d, length %d - wanted mode %d, length %d"),
229 :     TYPEOF(v), LENGTH(v), TYPEOF(vv), LENGTH(vv));
230 :     switch (TYPEOF(v)) {
231 :     case LGLSXP:
232 :     Memcpy(LOGICAL(vv), LOGICAL(v), LENGTH(vv));
233 :     break;
234 :     case INTSXP:
235 :     Memcpy(INTEGER(vv), INTEGER(v), LENGTH(vv));
236 :     break;
237 :     case REALSXP:
238 :     Memcpy(REAL(vv), REAL(v), LENGTH(vv));
239 :     break;
240 :     default:
241 :     error(_("invalid type for eval_check_store"));
242 :     }
243 :     UNPROTECT(1);
244 :     return vv;
245 : bates 415 }
246 :    
247 : bates 747 /**
248 : bates 1143 * Evaluate an expression in an environment, check that the length and
249 :     * mode are as expected and return the result.
250 : bates 747 *
251 : bates 1143 * @param fcn expression to evaluate
252 :     * @param rho environment in which to evaluate it
253 :     * @param mode desired mode
254 :     * @param len desired length
255 : bates 747 *
256 : bates 1143 * @return evaluated expression
257 : bates 747 */
258 : bates 1143 static SEXP
259 :     eval_check(SEXP fcn, SEXP rho, SEXPTYPE mode, int len) {
260 :     SEXP v = PROTECT(eval(fcn, rho));
261 :     if (TYPEOF(v) != mode || LENGTH(v) != len)
262 :     error(_("fcn produced mode %d, length %d - wanted mode %d, length %d"),
263 :     TYPEOF(v), LENGTH(v), mode, len);
264 :     UNPROTECT(1);
265 :     return v;
266 : bates 536 }
267 :    
268 : bates 1143 typedef struct glmer_struct
269 : bates 443 {
270 : bates 1143 SEXP cv; /* control values */
271 :     SEXP mer; /* mixed-effects representation */
272 :     SEXP rho; /* environment in which to evaluate the calls */
273 :     SEXP eta; /* linear predictor */
274 :     SEXP mu; /* mean vector */
275 :     SEXP linkinv; /* expression for inverse link evaluation */
276 :     SEXP mu_eta; /* expression for dmu/deta evaluation */
277 :     SEXP LMEopt; /* expression for LME optimization */
278 :     SEXP dev_resids; /* expression for deviance residuals */
279 :     SEXP var; /* expression for variance evaluation */
280 :     double *offset; /* offset for GLM */
281 :     double *wts; /* prior weights for GLM */
282 :     double *y; /* copy of response vector */
283 :     double *etaold; /* previous value of eta for evaluating */
284 :     int n; /* length of the response vector */
285 :     int p; /* length of fixed effects vector */
286 :     int nf; /* number of grouping factors */
287 :     int npar; /* total length of the parameter */
288 :     int niterEM; /* default number of ECME iterations */
289 :     int EMverbose; /* logical indicator */
290 :     int maxiter; /* maximum number of IRLS iterations */
291 :     double tol; /* convergence tolerance for IRLS iterations */
292 :     } glmer_struct, *GlmerStruct;
293 : bates 443
294 : bates 1143 /**
295 : bates 1162 * Calculate fitted values for the current fixed and random effects.
296 : bates 1143 *
297 :     * @param x pointer to an mer object
298 :     * @param initial initial values (i.e. an offset) or (double *) NULL
299 :     * @param val array to hold the fitted values
300 :     *
301 :     * @return pointer to a numeric array of fitted values
302 : bates 255 */
303 : bates 1143 static double *
304 : bates 1162 internal_mer_fitted(SEXP x, const double initial[], double val[])
305 : bates 255 {
306 : bates 1162 SEXP fixef = GET_SLOT(x, Matrix_fixefSym),
307 :     ranef = GET_SLOT(x, Matrix_ranefSym);
308 :     int ione = 1, n = LENGTH(GET_SLOT(x, Matrix_ySym)), p = LENGTH(fixef);
309 :     double *X = REAL(GET_SLOT(x, Matrix_XSym)), one[] = {1,0};
310 :     cholmod_sparse *Zt = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtSym));
311 :     cholmod_dense *chv = numeric_as_chm_dense(val, n),
312 :     *chb = as_cholmod_dense(ranef);
313 : bates 1143
314 :     mer_secondary(x);
315 :     if (initial) Memcpy(val, initial, n);
316 :     else AZERO(val, n);
317 : bates 1162 F77_CALL(dgemv)("N", &n, &p, one, X, &n, REAL(fixef), &ione, one, val, &ione);
318 :     if (!cholmod_sdmult(Zt, 1, one, one, chb, chv, &c))
319 :     error(_("Error return from sdmult"));
320 :     Free(chv); Free(chb); Free(Zt);
321 : bates 1143 return val;
322 : bates 255 }
323 : bates 1143
324 :     /**
325 :     * Extract the coefficients
326 :     *
327 :     * @param x pointer to an mer object
328 :     * @param ptyp parameter type to extract
329 :     * @param ans vector to hold the extracted values
330 :     *
331 :     * @return ans
332 : bates 255 */
333 : bates 813 static double *
334 : bates 1143 internal_mer_coef(SEXP x, int ptyp, double ans[])
335 : bates 255 {
336 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
337 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
338 : bates 813 i, nf = length(Omega), vind;
339 : bates 255
340 : bates 813 vind = 0; /* index in ans */
341 : bates 255 for (i = 0; i < nf; i++) {
342 :     int nci = nc[i], ncip1 = nci + 1;
343 :     if (nci == 1) {
344 : bates 1143 double dd = REAL(GET_SLOT(VECTOR_ELT(Omega, i), Matrix_xSym))[0];
345 : bates 813 ans[vind++] = ptyp ? ((ptyp == 1) ? log(dd) : 1./dd) : dd;
346 : bates 255 } else {
347 : bates 753 if (ptyp) { /* L log(D) L' factor of Omega[,,i] */
348 : bates 255 int j, k, ncisq = nci * nci;
349 : bates 1143 double *tmp = Memcpy(Calloc(ncisq, double),
350 :     REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)),
351 :     Matrix_xSym)), ncisq);
352 : bates 255 for (j = 0; j < nci; j++) {
353 :     double diagj = tmp[j * ncip1];
354 : bates 813 ans[vind++] = (ptyp == 1) ? (2. * log(diagj)) :
355 : bates 753 1./(diagj * diagj);
356 : bates 255 for (k = j + 1; k < nci; k++) {
357 :     tmp[j + k * nci] /= diagj;
358 :     }
359 :     }
360 :     for (j = 0; j < nci; j++) {
361 :     for (k = j + 1; k < nci; k++) {
362 : bates 813 ans[vind++] = tmp[j + k * nci];
363 : bates 255 }
364 :     }
365 :     Free(tmp);
366 :     } else { /* upper triangle of Omega[,,i] */
367 :     int j, k, odind = vind + nci;
368 : bates 1143 double *omgi = REAL(GET_SLOT(VECTOR_ELT(Omega, i), Matrix_xSym));
369 : bates 255
370 :     for (j = 0; j < nci; j++) {
371 : bates 813 ans[vind++] = omgi[j * ncip1];
372 : bates 255 for (k = j + 1; k < nci; k++) {
373 : bates 813 ans[odind++] = omgi[k*nci + j];
374 : bates 255 }
375 :     }
376 :     vind = odind;
377 :     }
378 :     }
379 :     }
380 : bates 813 return ans;
381 :     }
382 :    
383 : bates 1143 /**
384 :     * Set the coefficient vector and perform a factorization
385 : bates 813 *
386 : bates 1143 * @param x pointer to an mer object
387 :     * @param cc vector of coefficients to assign
388 :     * @param ptyp indicator of the parametrization being used
389 : bates 813 */
390 : bates 766 static
391 : bates 1143 void internal_mer_coefGets(SEXP x, const double cc[], int ptyp)
392 : bates 255 {
393 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
394 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
395 : bates 766 cind, i, nf = length(Omega);
396 : bates 255
397 :     cind = 0;
398 :     for (i = 0; i < nf; i++) {
399 : bates 1143 SEXP Omegai = VECTOR_ELT(Omega, i);
400 :     int j, k, nci = nc[i], ncisq = nc[i] * nc[i];
401 :     double *choli, *omgi = REAL(GET_SLOT(Omegai, Matrix_xSym));
402 :    
403 : bates 255 if (nci == 1) {
404 : bates 753 double dd = cc[cind++];
405 : bates 1143 *omgi = ptyp ? ((ptyp == 1) ? exp(dd) : 1./dd) : dd;
406 : bates 255 } else {
407 :     int odind = cind + nci, /* off-diagonal index */
408 : bates 1143 ncip1 = nci + 1;
409 :    
410 : bates 753 if (ptyp) {
411 : bates 1143 /* FIXME: Write this as an LDL decomposition */
412 : bates 753 double *tmp = Calloc(ncisq, double),
413 : bates 255 diagj, one = 1., zero = 0.;
414 :    
415 : bates 429 AZERO(omgi, ncisq);
416 : bates 255 for (j = 0; j < nci; j++) {
417 : bates 753 double dd = cc[cind++];
418 :     tmp[j * ncip1] = diagj =
419 :     (ptyp == 1) ? exp(dd/2.) : sqrt(1./dd);
420 : bates 255 for (k = j + 1; k < nci; k++) {
421 :     tmp[k*nci + j] = cc[odind++] * diagj;
422 :     }
423 :     }
424 :     F77_CALL(dsyrk)("U", "T", &nci, &nci, &one,
425 :     tmp, &nci, &zero, omgi, &nci);
426 :     Free(tmp);
427 :     } else {
428 :     for (j = 0; j < nci; j++) {
429 :     omgi[j * ncip1] = cc[cind++];
430 :     for (k = j + 1; k < nci; k++) {
431 :     omgi[k*nci + j] = cc[odind++];
432 :     }
433 :     }
434 :     }
435 :     cind = odind;
436 :     }
437 : bates 1143 choli = REAL(GET_SLOT(dpoMatrix_chol(Omegai), Matrix_xSym));
438 :     Memcpy(choli, omgi, ncisq);
439 :     F77_CALL(dpotrf)("U", &nci, choli, &nci, &j);
440 :     /* Yes, you really do need to do that decomposition.
441 :     The contents of choli before the decomposition are
442 :     from the previous value of Omegai. */
443 :     if (j)
444 :     error(_("Omega[[%d]] is not positive definite"), i + 1);
445 : bates 255 }
446 : bates 1143 flag_not_factored(x);
447 : bates 766 }
448 :    
449 : bates 1143 /**
450 :     * Evaluate current estimate of sigma from an mer object
451 :     *
452 :     * @param x pointer to an mer object
453 :     * @param REML indicator of whether to use REML.
454 :     * < 0 -> determine REML or ML from x@method
455 :     * == 0 -> use ML unconditionally
456 :     * > 0 -> use REML unconditionally
457 :     *
458 :     * @return
459 : bates 766 */
460 : bates 1143 static double
461 :     internal_mer_sigma(SEXP x, int REML)
462 : bates 766 {
463 : bates 1143 double *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym));
464 : bates 255
465 : bates 1143 if (REML < 0) /* get REML from x */
466 :     REML = !strcmp(CHAR(asChar(GET_SLOT(x,
467 :     Matrix_methodSym))),
468 :     "REML");
469 :     mer_factor(x);
470 :     return exp(dcmp[3]/2)/sqrt(dcmp[0] - (REML ? dcmp[1] : 0));
471 : bates 773 }
472 :    
473 : bates 1143 /**
474 :     * Update the derived quantities (ZtZ, ZtX, XtX, Zty, Xty
475 : bates 1156 * and dcmp[2] = y'y when Z, X, y, wts or wrkres has been changed.
476 : bates 1143 *
477 :     * @param x pointer to an mer object
478 :     * @param perm permutation from the Cholesky factor
479 : bates 255 */
480 :    
481 : bates 1143 static void
482 :     internal_mer_update_ZXy(SEXP x, int *perm)
483 : bates 255 {
484 : bates 1143 SEXP Xp = GET_SLOT(x, Matrix_XSym), ZtZ = GET_SLOT(x, Matrix_ZtZSym),
485 :     Ztyp = GET_SLOT(x, Matrix_ZtySym);
486 :     SEXP ZtZx = GET_SLOT(ZtZ, Matrix_xSym),
487 :     ZtZp = GET_SLOT(ZtZ, Matrix_pSym), ZtZi = GET_SLOT(ZtZ, Matrix_iSym);
488 :     int *dims = INTEGER(getAttrib(Xp, R_DimSymbol)), i, ione = 1, j;
489 :     int n = dims[0], nnz, p = dims[1], q = LENGTH(Ztyp);
490 :     cholmod_sparse *ts1, *ts2,
491 :     *Zt = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtSym));
492 : bates 1156 cholmod_sparse *Ztcp = cholmod_copy_sparse(Zt, &c);
493 :     int *Zp = (int*)Ztcp->p;
494 :     double *XtX = REAL(GET_SLOT(GET_SLOT(x, Matrix_XtXSym), Matrix_xSym)),
495 : bates 1143 *Xty = REAL(GET_SLOT(x, Matrix_XtySym)),
496 :     *ZtX = REAL(GET_SLOT(GET_SLOT(x, Matrix_ZtXSym), Matrix_xSym)),
497 :     *Zty = REAL(Ztyp),
498 : bates 1156 *wts = REAL(GET_SLOT(x, Matrix_wtsSym)),
499 : bates 1143 one[] = {1, 0}, zero[] = {0,0};
500 :     cholmod_dense *td1, *Xd = as_cholmod_dense(Xp),
501 : bates 1156 *wkrd = as_cholmod_dense(GET_SLOT(x, Matrix_wrkresSym));
502 :     cholmod_dense *Xcp = cholmod_copy_dense(Xd, &c),
503 :     *wkrcp = cholmod_copy_dense(wkrd, &c);
504 :     double *X = (double*)(Xcp->x), *Ztx = (double*)(Ztcp->x),
505 :     *wtres = (double*)(wkrcp->x);
506 :     /* Apply weights */
507 :     for (i = 0; i < n; i++)
508 :     wtres[i] = ((double*)(wkrd->x))[i] * wts[i];
509 :     for (j = 0; j < p; j++)
510 :     for (i = 0; i < n; i++)
511 :     X[i + j * n] = ((double*)(Xd->x))[i + j * n] * wts[i];
512 :     for (j = 0; j < n; j++)
513 :     for (i = Zp[j]; i < Zp[j + 1]; i++)
514 :     Ztx[i] = ((double*)(Zt->x))[i] * wts[j];
515 :     Free(Zt); Free(Xd); Free(wkrd);
516 :    
517 : bates 1143 /* y'y */
518 :     REAL(GET_SLOT(x, Matrix_devCompSym))[2] =
519 : bates 1156 F77_CALL(ddot)(&n, wtres, &ione, wtres, &ione);
520 : bates 1143 /* ZtZ */
521 : bates 1156 ts1 = cholmod_aat(Ztcp, (int *) NULL, (size_t) 0, 1/* mode */, &c);
522 : bates 1143 /* cholmod_aat returns stype == 0; copy to set stype == 1 */
523 :     ts2 = cholmod_copy(ts1, 1/* stype */, 1/* mode */, &c);
524 :     nnz = cholmod_nnz(ts2, &c);
525 :     if (((int)(ts2->ncol) + 1) != LENGTH(ZtZp))
526 :     error(_("Order of Z'Z has changed - was %d, now %d"),
527 :     LENGTH(ZtZp) - 1, (int)(ts2->ncol));
528 :     Memcpy(INTEGER(ZtZp), (int*)(ts2->p), LENGTH(ZtZp));
529 :     if (nnz != LENGTH(ZtZx))
530 :     error(_("Number of nonzeros in Z'Z has changed - was %d, now %d"),
531 :     LENGTH(ZtZx), nnz);
532 :     Memcpy(INTEGER(ZtZi), (int*)(ts2->i), nnz);
533 :     Memcpy(REAL(ZtZx), (double*)(ts2->x), nnz);
534 :     cholmod_free_sparse(&ts1, &c); cholmod_free_sparse(&ts2, &c);
535 :     /* PZ'X into ZtX */
536 :     td1 = cholmod_allocate_dense(q, p, q, CHOLMOD_REAL, &c);
537 : bates 1156 if (!cholmod_sdmult(Ztcp, 0, one, zero, Xcp, td1, &c))
538 : bates 1143 error(_("cholmod_sdmult failed"));
539 :     for (j = 0; j < p; j++) { /* apply the permutation to each column */
540 :     double *dcol = ZtX + j * q,
541 :     *scol = (double*)(td1->x) + j * q;
542 :     for (i = 0; i < q; i++) dcol[i] = scol[perm[i]];
543 : bates 255 }
544 : bates 1156 cholmod_free_dense(&td1, &c);
545 : bates 1143 /* PZ'y into Zty */
546 :     td1 = cholmod_allocate_dense(q, 1, q, CHOLMOD_REAL, &c);
547 : bates 1156 if (!cholmod_sdmult(Ztcp, 0, one, zero, wkrcp, td1, &c))
548 : bates 1143 error(_("cholmod_sdmult failed"));
549 :     for (i = 0; i < q; i++) Zty[i] = ((double *)(td1->x))[perm[i]];
550 : bates 1156 cholmod_free_dense(&td1, &c); cholmod_free_sparse(&Ztcp, &c);
551 : bates 1143 /* XtX and Xty */
552 :     AZERO(XtX, p * p);
553 :     F77_CALL(dsyrk)("U", "T", &p, &n, one, X, &n, zero, XtX, &p);
554 : bates 1156 F77_CALL(dgemv)("T", &n, &p, one, X, &n, wtres, &ione, zero, Xty, &ione);
555 :     cholmod_free_dense(&Xcp, &c); cholmod_free_dense(&wkrcp, &c);
556 : bates 1143 flag_not_factored(x);
557 : bates 255 }
558 :    
559 : bates 1143 static double chm_log_abs_det(cholmod_factor *F)
560 : bates 255 {
561 : bates 1143 double ans = 0;
562 : bates 255
563 : bates 1143 if (F->is_super) {
564 :     int i;
565 :     for (i = 0; i < F->nsuper; i++) {
566 :     int j, nrp1 = 1 + ((int *)(F->pi))[i + 1] - ((int *)(F->pi))[i],
567 :     nc = ((int *)(F->super))[i + 1] - ((int *)(F->super))[i];
568 :     double *x = (double *)(F->x) + ((int *)(F->px))[i];
569 : bates 422
570 : bates 1143 for (j = 0; j < nc; j++) ans += log(fabs(x[j * nrp1]));
571 : bates 255 }
572 : bates 1143 } else
573 :     error(_("code for simplicial factorization not yet written"));
574 :     return ans;
575 : bates 255 }
576 :    
577 : bates 1143 static double
578 :     Omega_log_det(SEXP Omega, int *nc, int *Gp)
579 : bates 255 {
580 : bates 1143 double ans = 0;
581 : bates 422 int i;
582 : bates 255
583 : bates 1143 for (i = 0; i < LENGTH(Omega); i++) {
584 :     int j, nci = nc[i], ncip1 = nc[i] + 1, nlev = (Gp[i + 1] - Gp[i])/nc[i];
585 :     double *omgi = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)),
586 :     Matrix_xSym));
587 : maechler 534
588 : bates 1143 for (j = 0; j < nci; j++) ans += 2. * nlev * log(fabs(omgi[j * ncip1]));
589 : bates 776 }
590 : bates 1143 return ans;
591 : bates 776 }
592 :    
593 : bates 255
594 : bates 761 /**
595 : bates 1143 * Inflate Z'Z to Z'Z+Omega and factor. Form RZX and rZy and update
596 :     * the status flags.
597 : bates 761 *
598 : bates 1143 * @param x pointer to an mer object.
599 : bates 761 */
600 : bates 1143 static void
601 :     internal_mer_Zfactor(SEXP x, cholmod_factor *L)
602 : bates 255 {
603 : bates 1143 SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),
604 :     rZyP = GET_SLOT(x, Matrix_rZySym);
605 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
606 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
607 :     nf = LENGTH(Omega), q = LENGTH(rZyP),
608 :     p = LENGTH(GET_SLOT(x, Matrix_rXySym));
609 :     cholmod_sparse *A, *Omg,
610 :     *zz = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtZSym));
611 :     cholmod_dense *ZtX = as_cholmod_dense(GET_SLOT(x, Matrix_ZtXSym)),
612 :     *Zty = as_cholmod_dense(GET_SLOT(x, Matrix_ZtySym)),
613 :     *rZy, *RZX;
614 :     int *omp, *nnz = Calloc(nf + 1, int), i,
615 :     *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
616 :     double *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym)), one[] = {1, 0};
617 :    
618 : bates 255
619 : bates 1143 dcmp[5] = Omega_log_det(Omega, nc, Gp); /* logDet(Omega) */
620 :     for (nnz[0] = 0, i = 0; i < nf; i++)
621 :     nnz[i + 1] = nnz[i] + ((Gp[i + 1] - Gp[i])*(nc[i] + 1))/2;
622 :     Omg = cholmod_allocate_sparse(zz->nrow, zz->ncol, (size_t) nnz[nf],
623 :     TRUE, TRUE, 1, CHOLMOD_REAL, &c);
624 :     omp = (int *) Omg->p;
625 : bates 255 for (i = 0; i < nf; i++) {
626 : bates 1143 int bb = Gp[i], j, jj, k, nci = nc[i];
627 :     int nlev = (Gp[i + 1] - bb)/nci;
628 :     double *Omgi = REAL(GET_SLOT(VECTOR_ELT(Omega, i), Matrix_xSym));
629 : bates 255
630 : bates 1143 for (j = 0; j < nlev; j++) { /* column of result */
631 :     int col0 = bb + j * nci; /* absolute column number */
632 : maechler 534
633 : bates 1143 for (jj = 0; jj < nci; jj++) { /* column of Omega_i */
634 :     int coljj = col0 + jj;
635 : bates 255
636 : bates 1143 omp[coljj + 1] = omp[coljj] + jj + 1;
637 :     for (k = 0; k <= jj; k++) { /* row of Omega_i */
638 :     int ind = omp[coljj];
639 :     ((int *)Omg->i)[ind + k] = col0 + k;
640 :     ((double *)Omg->x)[ind + k] = Omgi[jj * nci + k];
641 : bates 255 }
642 :     }
643 :     }
644 :     }
645 : bates 1143 Free(nnz);
646 :     A = cholmod_add(zz, Omg, one, one, TRUE, TRUE, &c);
647 :     Free(zz); cholmod_free_sparse(&Omg, &c);
648 :     if (!cholmod_factorize(A, L, &c))
649 :     error(_("rank_deficient Z'Z+Omega"));
650 :     cholmod_free_sparse(&A, &c);
651 :     dcmp[4] = 2 * chm_log_abs_det(L); /* 2 * logDet(L) */
652 :    
653 :     /* calculate and store RZX and rZy */
654 :     RZX = cholmod_solve(CHOLMOD_L, L, ZtX, &c); Free(ZtX);
655 :     rZy = cholmod_solve(CHOLMOD_L, L, Zty, &c); Free(Zty);
656 :     Memcpy(REAL(GET_SLOT(GET_SLOT(x, Matrix_RZXSym), Matrix_xSym)),
657 :     (double *) RZX->x, q * p);
658 :     cholmod_free_dense(&RZX, &c);
659 :     Memcpy(REAL(rZyP), (double *) rZy->x, q);
660 :     cholmod_free_dense(&rZy, &c);
661 :     /* signal that secondary slots are not valid */
662 :     status[0] = 1;
663 :     status[1] = status[2] = status[3] = 0;
664 : bates 255 }
665 :    
666 : bates 1143 /**
667 :     * Update the relative precision matrices by sampling from a Wishart
668 :     * distribution with scale factor determined by the current sample of
669 :     * random effects.
670 :     *
671 :     * @param Omega pointer to the list of relative precision matrices
672 :     * @param b current sample from the random effects
673 :     * @param sigma current value of sigma
674 :     * @param nf number of grouping factors
675 :     * @param nc number columns per grouping factor
676 :     * @param Gp integer vector pointing to the beginning of each outer
677 :     * group of columns
678 :     * @param vals vector in which to store values
679 :     * @param trans logical value indicating if variance components are
680 :     * stored in the transformed scale.
681 : bates 255 */
682 : bates 1143 static void
683 :     internal_Omega_update(SEXP Omega, const double b[], double sigma, int nf,
684 :     const int nc[], const int Gp[], double *vals, int trans)
685 : bates 255 {
686 : bates 1143 int i, j, k, info;
687 :     double one = 1, zero = 0;
688 : maechler 534
689 : bates 255 for (i = 0; i < nf; i++) {
690 : bates 824 int nci = nc[i];
691 : bates 1143 int nlev = (Gp[i + 1] - Gp[i])/nci, ncip1 = nci + 1,
692 :     ncisqr = nci * nci;
693 : bates 255 double
694 : bates 1143 *scal = Calloc(ncisqr, double), /* factor of scale matrix */
695 :     *tmp = Calloc(ncisqr, double),
696 :     *var = Calloc(ncisqr, double), /* simulated variance-covariance */
697 :     *wfac = Calloc(ncisqr, double); /* factor of Wishart variate */
698 : maechler 534
699 : bates 1143 /* generate and factor the scale matrix */
700 :     AZERO(scal, ncisqr);
701 :     F77_CALL(dsyrk)("U", "N", &nci, &nlev, &one, b + Gp[i], &nci,
702 :     &zero, scal, &nci);
703 :     F77_CALL(dpotrf)("U", &nci, scal, &nci, &info);
704 :     if (info)
705 :     error(_("Singular random effects varcov at level %d"), i + 1);
706 : bates 255
707 : bates 1143 /* generate a random factor from a standard Wishart distribution */
708 :     AZERO(wfac, ncisqr);
709 :     std_rWishart_factor((double) (nlev - nci + 1), nci, wfac);
710 : bates 255
711 : bates 1143 /* form the variance-covariance matrix and store elements */
712 :     Memcpy(tmp, scal, ncisqr);
713 :     F77_CALL(dtrsm)("L", "U", "T", "N", &nci, &nci,
714 :     &one, wfac, &nci, tmp, &nci);
715 :     F77_CALL(dsyrk)("U", "T", &nci, &nci, &one, tmp, &nci,
716 :     &zero, var, &nci);
717 :     for (j = 0; j < nci; j++) {
718 :     *vals++ = (trans ? log(var[j * ncip1]) : var[j * ncip1]);
719 : bates 775 }
720 : bates 1143 for (j = 1; j < nci; j++) {
721 :     for (k = 0; k < j; k++) {
722 :     *vals++ = (trans ? atanh(var[k + j * nci]/
723 :     sqrt(var[j * ncip1] * var[k * ncip1]))
724 :     : var[k + j * nci]);
725 : bates 766 }
726 :     }
727 : bates 1143 /* calculate and store the relative precision matrix */
728 :     Memcpy(tmp, wfac, ncisqr);
729 :     F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,
730 :     &sigma, scal, &nci, tmp, &nci);
731 :     F77_CALL(dsyrk)("U", "T", &nci, &nci, &one, tmp, &nci, &zero,
732 :     REAL(GET_SLOT(VECTOR_ELT(Omega, i), Matrix_xSym)),
733 :     &nci);
734 :     dpoMatrix_chol(VECTOR_ELT(Omega, i));
735 :     Free(scal); Free(tmp); Free(wfac); Free(var);
736 : bates 766 }
737 :     }
738 :    
739 : bates 818 /**
740 : bates 1143 * Evaluate new weights and working residuals.
741 : bates 818 *
742 : bates 776 * @param GS a GlmerStruct object
743 :     */
744 : bates 775 static void
745 : bates 1143 internal_glmer_reweight(GlmerStruct GS) {
746 : bates 1162 SEXP dmu_deta, var;
747 :     int i;
748 :     double *w = REAL(GET_SLOT(GS->mer, Matrix_wtsSym)),
749 :     *y = REAL(GET_SLOT(GS->mer, Matrix_ySym)),
750 :     *z = REAL(GET_SLOT(GS->mer, Matrix_wrkresSym));
751 : bates 1143
752 : bates 775 /* reweight mer */
753 :     eval_check_store(GS->linkinv, GS->rho, GS->mu);
754 :     dmu_deta = PROTECT(eval_check(GS->mu_eta, GS->rho,
755 :     REALSXP, GS->n));
756 :     var = PROTECT(eval_check(GS->var, GS->rho,
757 :     REALSXP, GS->n));
758 :     for (i = 0; i < GS->n; i++) {
759 : bates 1162 w[i] = GS->wts[i] *
760 : bates 775 REAL(dmu_deta)[i]/sqrt(REAL(var)[i]);
761 : bates 1162 z[i] = REAL(GS->eta)[i] - GS->offset[i] +
762 :     (y[i] - REAL(GS->mu)[i])/REAL(dmu_deta)[i];
763 : bates 775 }
764 :     UNPROTECT(2);
765 : bates 1143 mer_update_ZXy(GS->mer);
766 : bates 775 }
767 :    
768 : bates 776 /**
769 :     * Update eta, evaluate the convergence criterion, then copy eta to
770 :     * etaold
771 :     *
772 :     * @param GS a GlmerStruct object
773 :     * @param etaold previous values of the linear predictors
774 :     *
775 :     * @return convergence criterion
776 :     */
777 : bates 775 static double
778 : bates 795 conv_crit(double etaold[], double eta[], int n) {
779 : bates 782 double max_abs_eta = -1, max_abs_diff = -1;
780 : bates 775 int i;
781 :    
782 : bates 795 for (i = 0; i < n; i++) {
783 : bates 775 double abs_eta, abs_diff;
784 : bates 795
785 :     abs_eta = fabs(eta[i]);
786 : bates 782 if (abs_eta > max_abs_eta) max_abs_eta = abs_eta;
787 : bates 795 abs_diff = fabs(eta[i] - etaold[i]);
788 : bates 782 if (abs_diff > max_abs_diff) max_abs_diff = abs_diff;
789 : bates 795 etaold[i] = eta[i];
790 : bates 775 }
791 : bates 782 return max_abs_diff / (0.1 + max_abs_eta);
792 : bates 775 }
793 :    
794 : bates 776 /**
795 : bates 1143 * Update the ranef slot assuming that the fixef, rZy, RZX and L slots
796 :     * are up to date.
797 :     *
798 :     * @param x Pointer to an mer object
799 :     *
800 :     * @return
801 :     */
802 :     static double*
803 :     internal_mer_ranef(SEXP x)
804 :     {
805 :     SEXP ranef = GET_SLOT(x, Matrix_ranefSym);
806 :     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
807 :     if (!status[0]) {
808 :     error("Applying internal_mer_ranef without factoring");
809 :     return (double*)NULL; /* -Wall */
810 : bates 775 }
811 : bates 1143 if (!status[1]) {
812 :     SEXP fixef = GET_SLOT(x, Matrix_fixefSym),
813 :     ranef = GET_SLOT(x, Matrix_ranefSym);
814 :     int ione = 1, p = LENGTH(fixef), q = LENGTH(ranef);
815 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
816 :     cholmod_dense *td1, *td2,
817 :     *chranef = as_cholmod_dense(ranef);
818 :     double *RZX = REAL(GET_SLOT(GET_SLOT(x, Matrix_RZXSym), Matrix_xSym)),
819 :     m1[] = {-1,0}, one[] = {1,0};
820 :    
821 :     Memcpy(REAL(ranef), REAL(GET_SLOT(x, Matrix_rZySym)), q);
822 :     F77_CALL(dgemv)("N", &q, &p, m1, RZX, &q, REAL(fixef), &ione,
823 :     one, REAL(ranef), &ione);
824 :     td1 = cholmod_solve(CHOLMOD_Lt, L, chranef, &c);
825 :     td2 = cholmod_solve(CHOLMOD_Pt, L, td1, &c);
826 :     Free(chranef); cholmod_free_dense(&td1, &c);
827 :     Memcpy(REAL(ranef), (double *)(td2->x), q);
828 :     cholmod_free_dense(&td2, &c);
829 :     status[1] = 1;
830 :     status[2] = status[3] = 0;
831 :     Free(L);
832 :     }
833 :     return REAL(ranef);
834 :     }
835 : bates 775
836 : bates 1143 /**
837 :     * Update the fixef slot on a factored mer object.
838 :     *
839 :     * @param x Pointer to an mer object
840 :     *
841 :     * @return fixed effects vector
842 :     */
843 :     static double*
844 :     internal_mer_fixef(SEXP x)
845 :     {
846 :     SEXP fixef = GET_SLOT(x, Matrix_fixefSym);
847 :     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
848 :     if (!status[0]) {
849 :     error("Applying internal_mer_fixef without factoring");
850 :     return (double*)NULL; /* -Wall */
851 :     }
852 :     if (!status[1]) {
853 :     int ione = 1, p = LENGTH(fixef);
854 :     Memcpy(REAL(fixef), REAL(GET_SLOT(x, Matrix_rXySym)), p);
855 :     F77_CALL(dtrsv)("U", "N", "N", &p,
856 :     REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym),
857 :     Matrix_xSym)),
858 :     &p, REAL(fixef), &ione);
859 :     }
860 :     return REAL(fixef);
861 : bates 775 }
862 : bates 1143
863 : bates 775 /**
864 : bates 1162 * Iterate to determine the conditional modes of the random effects.
865 : bates 787 *
866 :     * @param GS a GlmerStruct object
867 : bates 1143 * @param fixed vector of fixed effects
868 :     * @param varc vector of parameters for the variance-covariance
869 : bates 818 *
870 :     * @return An indicator of whether the iterations converged
871 : bates 787 */
872 : bates 804 static int
873 : bates 787 internal_bhat(GlmerStruct GS, const double fixed[], const double varc[])
874 :     {
875 : bates 1143 SEXP fixef = GET_SLOT(GS->mer, Matrix_fixefSym);
876 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(GS->mer, Matrix_LSym));
877 :     int i;
878 :     double crit = GS->tol + 1;
879 :    
880 :     if (varc) /* skip this step if varc == (double*) NULL */
881 :     internal_mer_coefGets(GS->mer, varc, 2);
882 :     Memcpy(REAL(fixef), fixed, LENGTH(fixef));
883 :     internal_mer_Zfactor(GS->mer, L);
884 :     internal_mer_ranef(GS->mer);
885 : bates 1162 internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
886 : bates 1143 Memcpy(GS->etaold, REAL(GS->eta), GS->n);
887 :    
888 :     for (i = 0; i < GS->maxiter && crit > GS->tol; i++) {
889 :     internal_glmer_reweight(GS);
890 :     internal_mer_Zfactor(GS->mer, L);
891 :     internal_mer_ranef(GS->mer);
892 : bates 1162 internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
893 : bates 1143 crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);
894 :     }
895 :     Free(L);
896 :     return (crit > GS->tol) ? 0 : i;
897 :     }
898 : bates 787
899 : bates 1143 /**
900 :     * Print the verbose output in the ECME iterations
901 :     *
902 :     * @param x pointer to an ssclme object
903 :     * @param iter iteration number
904 :     * @param REML non-zero for REML, zero for ML
905 :     */
906 :     static void
907 :     EMsteps_verbose_print(SEXP x, int iter, int REML)
908 :     {
909 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),
910 :     gradComp = GET_SLOT(x, Matrix_gradCompSym);
911 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
912 :     i, ifour = 4, ii, ione = 1, jj,
913 :     nf = LENGTH(Omega);
914 :     double
915 :     cc[] = {-1, 1, 1, REML ? 1 : 0},
916 :     *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),
917 :     one = 1., zero = 0.;
918 : bates 813
919 : bates 1143 if (iter == 0) Rprintf(" EM iterations\n");
920 :     Rprintf("%3d %.3f", iter, dev[REML ? 1 : 0]);
921 :     for (i = 0; i < nf; i++) {
922 :     int nci = nc[i], ncip1 = nci + 1, ncisqr = nci * nci;
923 :     double
924 :     *Omgi = REAL(GET_SLOT(VECTOR_ELT(Omega, i), Matrix_xSym)),
925 :     *Grad = Calloc(ncisqr, double);
926 : bates 787
927 : bates 1143 /* diagonals */
928 :     Rprintf(" (%#8g", Omgi[0]);
929 :     for (jj = 1; jj < nci; jj++) {
930 :     Rprintf(" %#8g", Omgi[jj * ncip1]);
931 :     }
932 :     for (jj = 1; jj < nci; jj++) /* offdiagonals */
933 :     for (ii = 0; ii < jj; ii++)
934 :     Rprintf(" %#8g", Omgi[ii + jj * nci]);
935 :     /* Evaluate and print the gradient */
936 :     F77_CALL(dgemv)("N", &ncisqr, &ifour, &one,
937 :     REAL(VECTOR_ELT(gradComp, i)), &ncisqr,
938 :     cc, &ione, &zero, Grad, &ione);
939 :     Rprintf(":%#8.3g", Grad[0]);
940 :    
941 :     for (jj = 1; jj < nci; jj++) { /* diagonals */
942 :     Rprintf(" %#8.3g", Grad[jj * ncip1]);
943 :     }
944 :     for (jj = 1; jj < nci; jj++) /* offdiagonals */
945 :     for (ii = 0; ii < jj; ii++)
946 :     Rprintf(" %#8.3g", Grad[ii + jj * nci]);
947 :     Rprintf(")");
948 :     Free(Grad);
949 : bates 787 }
950 : bates 1143 Rprintf("\n");
951 : bates 787 }
952 :    
953 :     /**
954 : bates 1143 * Perform a number of ECME steps
955 :     *
956 :     * @param x pointer to an mer object
957 :     * @param nEM number of iteration to perform
958 :     * @param verb indicator of verbose output
959 :     */
960 :     static void
961 :     internal_ECMEsteps(SEXP x, int nEM, int verb)
962 :     {
963 :     SEXP Omega = GET_SLOT(x, Matrix_OmegaSym),
964 :     gradComp = GET_SLOT(x, Matrix_gradCompSym);
965 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
966 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
967 :     REML = !strcmp(CHAR(asChar(GET_SLOT(x, Matrix_methodSym))),
968 :     "REML"),
969 :     i, ifour = 4, info, ione = 1, iter,
970 :     nf = LENGTH(Omega);
971 :     double
972 :     cc[] = {0, 1, 1, REML ? 1 : 0},
973 :     zero = 0.0;
974 :    
975 :     mer_gradComp(x);
976 :     if (verb)
977 :     EMsteps_verbose_print(x, 0, REML);
978 :     for (iter = 0; iter < nEM; iter++) {
979 :     for (i = 0; i < nf; i++) {
980 :     SEXP Omgi = VECTOR_ELT(Omega, i);
981 :     int nci = nc[i], ncisqr = nci * nci,
982 :     nlev = (Gp[i + 1] - Gp[i])/nci;
983 :     double *Omgx = REAL(GET_SLOT(Omgi, Matrix_xSym)),
984 :     mult = 1./((double) nlev);
985 :    
986 :     F77_CALL(dgemm)("N", "N", &ncisqr, &ione, &ifour, &mult,
987 :     REAL(VECTOR_ELT(gradComp, i)), &ncisqr,
988 :     cc, &ifour, &zero, Omgx, &ncisqr);
989 :     F77_CALL(dpotrf)("U", &nci, Omgx, &nci, &info);
990 :     if (info)
991 :     error(_("DPOTRF in ECME update gave code %d"),
992 :     info);
993 :     F77_CALL(dpotri)("U", &nci, Omgx, &nci, &info);
994 :     if (info)
995 :     error(_("Matrix inverse in ECME update gave code %d"), info);
996 :     SET_SLOT(Omgi, Matrix_factorSym, allocVector(VECSXP, 0));
997 :     }
998 :     flag_not_factored(x);
999 :     mer_gradComp(x);
1000 :     if (verb)
1001 :     EMsteps_verbose_print(x, iter + 1, REML);
1002 :     }
1003 :     }
1004 :    
1005 :     /**
1006 :     * Evaluate the conditional deviance for the stored random effects.
1007 :     *
1008 :     * @param GS Pointer to a GlmerStruct
1009 :     * @param fixed value of the fixed effects
1010 :     *
1011 :     * @return conditional deviance
1012 :     */
1013 :     static double
1014 :     random_effects_deviance(GlmerStruct GS)
1015 :     {
1016 :     SEXP devs;
1017 :     int i;
1018 :     double ans;
1019 :    
1020 : bates 1162 internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
1021 : bates 1143 eval_check_store(GS->linkinv, GS->rho, GS->mu);
1022 :     devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));
1023 :     for (i = 0, ans = 0; i < GS->n; i++) ans += REAL(devs)[i];
1024 :     UNPROTECT(1);
1025 :     return ans;
1026 :     }
1027 :    
1028 :     /* Externally accessible functions */
1029 :    
1030 : bates 786 /**
1031 : bates 1143 * Simulate a sample of random matrices from a Wishart distribution
1032 :     *
1033 :     * @param ns Number of samples to generate
1034 :     * @param dfp Degrees of freedom
1035 :     * @param scal Positive-definite scale matrix
1036 :     *
1037 :     * @return
1038 :     */
1039 :     SEXP
1040 :     Matrix_rWishart(SEXP ns, SEXP dfp, SEXP scal)
1041 :     {
1042 :     SEXP ans;
1043 :     int *dims = INTEGER(getAttrib(scal, R_DimSymbol)), j,
1044 :     n = asInteger(ns), psqr;
1045 :     double *scCp, *tmp, df = asReal(dfp), one = 1, zero = 0;
1046 :    
1047 :     if (!isMatrix(scal) || !isReal(scal) || dims[0] != dims[1])
1048 :     error("scal must be a square, real matrix");
1049 :     if (n <= 0) n = 1;
1050 :     psqr = dims[0] * dims[0];
1051 :     tmp = Calloc(psqr, double);
1052 :     AZERO(tmp, psqr);
1053 :     scCp = Memcpy(Calloc(psqr, double), REAL(scal), psqr);
1054 :     F77_CALL(dpotrf)("U", &(dims[0]), scCp, &(dims[0]), &j);
1055 :     if (j)
1056 :     error("scal matrix is not positive-definite");
1057 :     PROTECT(ans = alloc3Darray(REALSXP, dims[0], dims[0], n));
1058 :    
1059 :     GetRNGstate();
1060 :     for (j = 0; j < n; j++) {
1061 :     double *ansj = REAL(ans) + j * psqr;
1062 :     std_rWishart_factor(df, dims[0], tmp);
1063 :     F77_CALL(dtrmm)("R", "U", "N", "N", dims, dims,
1064 :     &one, scCp, dims, tmp, dims);
1065 :     F77_CALL(dsyrk)("U", "T", &(dims[1]), &(dims[1]),
1066 :     &one, tmp, &(dims[1]),
1067 :     &zero, ansj, &(dims[1]));
1068 :     internal_symmetrize(ansj, dims[0]);
1069 :     }
1070 :    
1071 :     PutRNGstate();
1072 :     Free(scCp); Free(tmp);
1073 :     UNPROTECT(1);
1074 :     return ans;
1075 :     }
1076 :    
1077 :     /**
1078 :     * Perform the PQL optimization
1079 :     *
1080 :     * @param GSp pointer to a GlmerStruct object
1081 :     *
1082 :     * @return R_NilValue
1083 :     */
1084 :     SEXP glmer_PQL(SEXP GSp)
1085 :     {
1086 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
1087 :     int i; double crit;
1088 :    
1089 :     Memcpy(GS->etaold, REAL(GS->eta), GS->n);
1090 :     for (i = 0, crit = GS->tol + 1;
1091 :     i < GS->maxiter && crit > GS->tol; i++) {
1092 :     internal_glmer_reweight(GS);
1093 :     if (!i) mer_initial(GS->mer); /* initialize first fit */
1094 :     internal_ECMEsteps(GS->mer, i ? 2 : GS->niterEM,
1095 :     GS->EMverbose);
1096 :     eval(GS->LMEopt, GS->rho);
1097 : bates 1162 internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
1098 : bates 1143 crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);
1099 :     }
1100 :     if (crit > GS->tol)
1101 :     warning(_("IRLS iterations for PQL did not converge"));
1102 :    
1103 :     return R_NilValue;
1104 :     }
1105 :    
1106 :     /**
1107 :     * Compute the Laplace approximation to the deviance.
1108 :     *
1109 :     * @param pars pointer to a numeric vector of parameters
1110 :     * @param GSp pointer to a GlmerStruct object
1111 :     *
1112 :     * @return the Laplace approximation to the deviance
1113 :     */
1114 :     SEXP glmer_devLaplace(SEXP pars, SEXP GSp)
1115 :     {
1116 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
1117 :     SEXP Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);
1118 :     int *Gp = INTEGER(GET_SLOT(GS->mer, Matrix_GpSym)),
1119 :     *nc = INTEGER(GET_SLOT(GS->mer, Matrix_ncSym));
1120 :     double *bhat = REAL(GET_SLOT(GS->mer, Matrix_ranefSym)),
1121 :     *dcmp = REAL(GET_SLOT(GS->mer, Matrix_devCompSym)),
1122 :     *dev = REAL(GET_SLOT(GS->mer, Matrix_devianceSym));
1123 :    
1124 :     if (!isReal(pars) || LENGTH(pars) != GS->npar)
1125 :     error(_("`%s' must be a numeric vector of length %d"),
1126 :     "pars", GS->npar);
1127 :     if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))
1128 :     return ScalarReal(DBL_MAX);
1129 :     dev[0] = dcmp[4] - dcmp[5] + random_effects_deviance(GS) +
1130 :     b_quadratic(bhat, Omega, Gp, nc);
1131 :     dev[1] = NA_REAL;
1132 :     return ScalarReal(dev[0]);
1133 :     }
1134 :    
1135 :     /**
1136 :     * Release the storage for a GlmerStruct
1137 : bates 805 *
1138 : bates 1143 * @param GSp External pointer to a GlmerStruct
1139 : bates 805 *
1140 : bates 1143 * @return R_NilValue
1141 : bates 805 */
1142 : bates 1143 SEXP glmer_finalize(SEXP GSp) {
1143 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
1144 :    
1145 : bates 1162 Free(GS->offset); Free(GS->wts); Free(GS->etaold);
1146 : bates 1143 Free(GS);
1147 :     return R_NilValue;
1148 : bates 805 }
1149 :    
1150 : bates 1143 /**
1151 :     * Return an external pointer object to a GlmerStruct created in
1152 :     * environment rho
1153 :     *
1154 :     * @param rho An environment
1155 :     *
1156 :     * @return An external pointer to a GlmerStruct
1157 :     */
1158 :     SEXP glmer_init(SEXP rho) {
1159 :     GlmerStruct GS;
1160 :     SEXP tmp, y, Ztx;
1161 :    
1162 :    
1163 :     GS = (GlmerStruct) Calloc(1, glmer_struct);
1164 :     if (!isEnvironment(rho))
1165 :     error(_("`rho' must be an environment"));
1166 :     GS->rho = rho;
1167 :     GS->mer = find_and_check(rho, install("mer"), VECSXP, 0);
1168 :     y = GET_SLOT(GS->mer, Matrix_ySym);
1169 :     GS->n = LENGTH(y);
1170 :     GS->p = LENGTH(GET_SLOT(GS->mer, Matrix_rXySym));
1171 :     GS->y = Memcpy(Calloc(GS->n, double), REAL(y), GS->n);
1172 :     Ztx = GET_SLOT(GET_SLOT(GS->mer, Matrix_ZtSym), Matrix_xSym);
1173 :     GS->mu = find_and_check(rho, install("mu"), REALSXP, GS->n);
1174 :     tmp = find_and_check(rho, install("offset"), REALSXP, GS->n);
1175 :     GS->offset = Memcpy(Calloc(GS->n, double), REAL(tmp), GS->n);
1176 :     tmp = find_and_check(rho, install("wts"), REALSXP, GS->n);
1177 :     GS->wts = Memcpy(Calloc(GS->n, double), REAL(tmp), GS->n);
1178 :     GS->etaold = Calloc(GS->n, double);
1179 :     GS->cv = find_and_check(rho, install("cv"), VECSXP, 0);
1180 :     GS->niterEM = asInteger(Matrix_getElement(GS->cv, "niterEM"));
1181 :     GS->EMverbose = asLogical(Matrix_getElement(GS->cv, "EMverbose"));
1182 :     GS->tol = asReal(Matrix_getElement(GS->cv, "tolerance"));
1183 :     GS->maxiter = asInteger(Matrix_getElement(GS->cv, "maxIter"));
1184 :     GS->nf = LENGTH(GET_SLOT(GS->mer, Matrix_flistSym));
1185 :     GS->npar = GS->p +
1186 :     coef_length(GS->nf, INTEGER(GET_SLOT(GS->mer, Matrix_ncSym)));
1187 :     GS->eta = find_and_check(rho, install("eta"), REALSXP, GS->n);
1188 : bates 805
1189 : bates 1143 GS->linkinv = find_and_check(rho, install("linkinv"),
1190 :     LANGSXP, 0);
1191 :     GS->mu_eta = find_and_check(rho, install("mu.eta"),
1192 :     LANGSXP, 0);
1193 :     GS->var = find_and_check(rho, install("variance"),
1194 :     LANGSXP, 0);
1195 :     GS->LMEopt = find_and_check(rho, install("doLMEopt"),
1196 :     LANGSXP, 0);
1197 :     GS->dev_resids = find_and_check(rho, install("dev.resids"),
1198 :     LANGSXP, 0);
1199 :     return R_MakeExternalPtr(GS, R_NilValue, GS->mer);
1200 : bates 781 }
1201 :    
1202 : bates 1143 /**
1203 :     * Perform ECME steps for the REML or ML criterion.
1204 :     *
1205 :     * @param x pointer to an mer object
1206 :     * @param nsteps pointer to an integer scalar - the number of ECME
1207 :     * steps to perform
1208 :     * @param Verbp pointer to a logical scalar indicating verbose output
1209 :     *
1210 :     * @return R_NilValue
1211 : bates 802 */
1212 : bates 1143 SEXP mer_ECMEsteps(SEXP x, SEXP nsteps, SEXP Verbp)
1213 : bates 802 {
1214 : bates 1143 int nstp = asInteger(nsteps);
1215 :     if (nstp > 0) internal_ECMEsteps(x, nstp, asLogical(Verbp));
1216 :     return R_NilValue;
1217 :     }
1218 : bates 802
1219 : bates 1143 /**
1220 :     * Fill in five symmetric matrices, providing the information to
1221 :     * generate the Hessian.
1222 :     *
1223 :     * @param x pointer to an mer object
1224 :     *
1225 :     * @return an array consisting of five symmetric faces
1226 : bates 804 */
1227 : bates 1143 SEXP mer_Hessian(SEXP x)
1228 : bates 802 {
1229 : bates 1143 SEXP
1230 :     D = GET_SLOT(x, Matrix_DSym),
1231 :     Omega = GET_SLOT(x, Matrix_OmegaSym),
1232 :     RZXP = GET_SLOT(x, Matrix_RZXSym),
1233 :     gradComp = GET_SLOT(x, Matrix_gradCompSym),
1234 :     val;
1235 :     int *dRZX = INTEGER(getAttrib(RZXP, R_DimSymbol)),
1236 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1237 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1238 :     Q, Qsqr, RZXpos, facepos,
1239 :     i, ione = 1, j, nf = length(Omega), p = dRZX[1] - 1, pos;
1240 :     double
1241 :     *RZX = REAL(RZXP),
1242 :     *b = REAL(RZXP) + dRZX[0] * p,
1243 :     *bbface, /* vec of second faces of gradComp elts */
1244 :     one = 1.,
1245 :     zero = 0.;
1246 : bates 802
1247 : bates 1143 mer_gradComp(x);
1248 :     Q = 0; /* number of rows and columns in the result */
1249 :     for (i = 0; i < nf; i++) Q += nc[i] * nc[i];
1250 :     Qsqr = Q * Q;
1251 :     bbface = Calloc(Q, double);
1252 :     val = PROTECT(alloc3Darray(REALSXP, Q, Q, 5));
1253 :     AZERO(REAL(val), Qsqr * 5);
1254 :    
1255 :     pos = 0;
1256 :     for (i = 0; i < nf; i++) {
1257 :     int nci = nc[i];
1258 :     int ncisqr = nci * nci;
1259 :     double *fDi = REAL(VECTOR_ELT(gradComp, i)),
1260 :     mult = 1./((double)(Gp[i + 1] - Gp[i])/nci);
1261 :    
1262 :     Memcpy(bbface + pos, fDi + ncisqr, ncisqr);
1263 :     /* outer product of the third face of gradComp on the diagonal
1264 :     * of the third face of val */
1265 :     F77_CALL(dsyr)("U", &ncisqr, &mult, fDi + 2 * ncisqr, &ione,
1266 :     REAL(val) + 2 * Qsqr + pos * Q, &Q);
1267 :     pos += ncisqr;
1268 : bates 784 }
1269 : bates 1143 /* fifth face is outer product of bbface */
1270 :     F77_CALL(dsyr)("U", &Q, &one, bbface, &ione, REAL(val) + 4 * Qsqr, &Q);
1271 :     /* fourth face from \bb\trans\der\vb\der\bb */
1272 :     AZERO(REAL(val) + 3 * Qsqr, Qsqr); /* zero accumulator */
1273 :     RZXpos = 0;
1274 :     facepos = 0;
1275 :     for (i = 0; i < nf; i++) {
1276 :     int ii, jj, nci = nc[i];
1277 :     int ncisqr = nci * nci, nctp = nci * p,
1278 :     nlev = (Gp[i + 1] - Gp[i])/nci;
1279 :     int maxpq = (p > nci) ? p : nci;
1280 :     double
1281 :     *Di = REAL(VECTOR_ELT(D, i)),
1282 :     *arr = Calloc(ncisqr * maxpq, double), /* tmp 3Darray */
1283 :     *face = REAL(val) + 3 * Qsqr,
1284 :     *mat = Calloc(nci * maxpq, double); /* tmp matrix */
1285 :    
1286 :     for (j = 0; j < nlev; j++) {
1287 :     F77_CALL(dgemm)("T", "T", &p, &nci, &nci,
1288 :     &one, RZX + j * nci, dRZX, Di + j * ncisqr, &nci,
1289 :     &zero, mat, &p);
1290 :     F77_CALL(dgemm)("N", "N", &nctp, &nci, &ione,
1291 :     &one, mat, &nctp, b + j * nci, &ione,
1292 :     &zero, arr, &nctp);
1293 :     F77_CALL(dsyrk)("U", "T", &ncisqr, &p, &one, arr, &p,
1294 :     &one, face + facepos, &Q);
1295 :     /* Add the D_{i,j}^{-T/2} term */
1296 :     Memcpy(mat, Di + j * ncisqr, ncisqr);
1297 :     for (jj = 1; jj < nci; jj++) { /* transpose mat */
1298 :     for (ii = 0; ii < jj; ii++) {
1299 :     mat[jj + ii * nci] = mat[ii + jj * nci];
1300 :     mat[ii + jj * nci] = 0.;
1301 :     }
1302 :     }
1303 :     F77_CALL(dgemm)("N", "N", &ncisqr, &nci, &ione,
1304 :     &one, mat, &ncisqr, b + j * nci, &ione,
1305 :     &zero, arr, &ncisqr);
1306 :     /* FIXME: Next call could be dsyr (it's rank one). */
1307 :     F77_CALL(dsyrk)("U", "T", &ncisqr, &nci, &one, arr, &nci,
1308 :     &one, face + facepos, &Q);
1309 :    
1310 : bates 781 }
1311 : bates 1143 RZXpos += nci * nlev;
1312 :     facepos += ncisqr;
1313 :     Free(arr); Free(mat);
1314 : bates 781 }
1315 : bates 1143 UNPROTECT(2);
1316 :     Free(bbface);
1317 :     return val;
1318 : bates 781 }
1319 : bates 818
1320 : bates 814 /**
1321 : bates 1143 * Generate a Markov-Chain Monte Carlo sample from a fitted
1322 :     * linear mixed model.
1323 : bates 814 *
1324 : bates 1143 * @param mer pointer to an mer object
1325 :     * @param savebp pointer to a logical scalar indicating if the
1326 :     * random-effects should be saved
1327 :     * @param nsampp pointer to an integer scalar of the number of samples
1328 :     * to generate
1329 :     * @param transp pointer to an logical scalar indicating if the
1330 :     * variance components should be transformed.
1331 : bates 814 *
1332 : bates 1143 * @return a matrix
1333 : bates 814 */
1334 : bates 1143 SEXP
1335 :     mer_MCMCsamp(SEXP x, SEXP savebp, SEXP nsampp, SEXP transp)
1336 : bates 814 {
1337 : bates 1143 SEXP ans, Omega = GET_SLOT(x, Matrix_OmegaSym),
1338 :     Omegacp = PROTECT(duplicate(Omega));
1339 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
1340 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1341 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1342 :     REML = !strcmp(CHAR(asChar(GET_SLOT(x, Matrix_methodSym))), "REML"),
1343 :     i, ione = 1, j, n = LENGTH(GET_SLOT(x, Matrix_ySym)),
1344 :     nf = LENGTH(Omega), nsamp = asInteger(nsampp),
1345 :     p = LENGTH(GET_SLOT(x, Matrix_rXySym)),
1346 :     q = LENGTH(GET_SLOT(x, Matrix_rZySym)),
1347 :     saveb = asLogical(savebp),
1348 :     trans = asLogical(transp);
1349 :     double
1350 :     *RXX = REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym), Matrix_xSym)),
1351 :     *RZX = REAL(GET_SLOT(GET_SLOT(x, Matrix_RZXSym), Matrix_xSym)),
1352 :     *bhat = REAL(GET_SLOT(x, Matrix_ranefSym)),
1353 :     *betahat = REAL(GET_SLOT(x, Matrix_fixefSym)),
1354 :     *bnew = Calloc(q, double), *betanew = Calloc(p, double),
1355 :     *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym)),
1356 :     df = n - (REML ? p : 0), m1[] = {-1,0}, one[] = {1,0};
1357 :     int nrbase = p + 1 + coef_length(nf, nc); /* rows always included */
1358 :     int nrtot = nrbase + (saveb ? q : 0);
1359 :     cholmod_dense *chb, *chbnew = numeric_as_chm_dense(bnew, q);
1360 :    
1361 :     if (nsamp <= 0) nsamp = 1;
1362 :     ans = PROTECT(allocMatrix(REALSXP, nrtot, nsamp));
1363 : bates 814 GetRNGstate();
1364 : bates 1143 for (i = 0; i < nsamp; i++) {
1365 :     double *col = REAL(ans) + i * nrtot, sigma;
1366 :     /* factor x and get secondary values */
1367 :     mer_factor(x);
1368 :     mer_secondary(x);
1369 :     /* simulate and store new value of sigma */
1370 :     sigma = exp(dcmp[3]/2)/sqrt(rchisq(df));
1371 :     col[p] = (trans ? 2 * log(sigma) : sigma * sigma);
1372 :     /* simulate scaled, independent normals */
1373 :     for (j = 0; j < p; j++) betanew[j] = sigma * norm_rand();
1374 :     for (j = 0; j < q; j++) bnew[j] = sigma * norm_rand();
1375 :     /* betanew := RXX^{-1} %*% betanew */
1376 :     F77_CALL(dtrsv)("U", "N", "N", &p, RXX, &p, betanew, &ione);
1377 :     /* bnew := bnew - RZX %*% betanew */
1378 :     F77_CALL(dgemv)("N", &q, &p, m1, RZX, &q, betanew, &ione,
1379 :     one, bnew, &ione);
1380 :     /* chb := L^{-T} %*% bnew */
1381 :     chb = cholmod_solve(CHOLMOD_Lt, L, chbnew, &c);
1382 :     /* Copy chb to bnew and free chb */
1383 :     Memcpy(bnew, (double *)(chb->x), q);
1384 :     cholmod_free_dense(&chb, &c);
1385 :     /* Add conditional modes and store beta */
1386 :     for (j = 0; j < p; j++) {
1387 :     col[j] = (betanew[j] += betahat[j]);
1388 :     }
1389 :     /* Add conditional modes and
1390 :     * optionally store b */
1391 :     for (j = 0; j < q; j++) {
1392 :     bnew[j] += bhat[j];
1393 :     if (saveb) col[nrbase + j] = bnew[j];
1394 :     }
1395 :     /* Update and store variance-covariance of the random effects */
1396 :     internal_Omega_update(Omega, bnew, sigma, nf, nc, Gp, col + p + 1, trans);
1397 :     }
1398 : bates 814 PutRNGstate();
1399 : bates 1143 Free(betanew); Free(bnew); Free(chbnew);
1400 :     /* Restore original Omega */
1401 :     SET_SLOT(x, Matrix_OmegaSym, Omegacp);
1402 :     mer_factor(x);
1403 : bates 812
1404 : bates 1143 Free(L);
1405 :     UNPROTECT(2);
1406 : bates 812 return ans;
1407 :     }
1408 :    
1409 : bates 1143 /**
1410 :     * Create an mer object from a list of grouping factors and a list of model
1411 :     * matrices.
1412 :     *
1413 :     * @param fl named list of grouping factors
1414 :     * @param Ztl list of transposes of model matrices
1415 :     * @param Xp model matrix for the fixed effects
1416 :     * @param yp response vector
1417 :     * @param method character vector describing the estimation method
1418 :     *
1419 :     * @return pointer to an mer object
1420 : bates 818 */
1421 : bates 1143 SEXP mer_create(SEXP fl, SEXP ZZt, SEXP Xp, SEXP yp, SEXP method,
1422 :     SEXP ncp, SEXP cnames, SEXP useS, SEXP call,
1423 :     SEXP family)
1424 : bates 812 {
1425 : bates 1143 SEXP Omega, bVar, gradComp, fnms = getAttrib(fl, R_NamesSymbol),
1426 :     val = PROTECT(NEW_OBJECT(MAKE_CLASS("mer"))), xnms;
1427 :     cholmod_sparse *ts1, *ts2, *Zt;
1428 :     cholmod_factor *F;
1429 :     int *nc = INTEGER(ncp), *Gp, *xdims, i,
1430 :     nf = LENGTH(fl), nobs = LENGTH(yp), p, q;
1431 :     char *devCmpnms[] = {"n", "p", "yty", "logryy2", "logDetL2",
1432 :     "logDetOmega", "logDetRXX", ""};
1433 :     char *devnms[] = {"ML", "REML", ""};
1434 :     char *statusnms[] =
1435 :     {"factored", "secondary", "gradComp", "Hesscomp", ""};
1436 : bates 1156 double *dcmp, *wts, *wrkres, *y;
1437 : bates 1143 /* Check arguments to be duplicated */
1438 :     if (!isReal(yp)) error(_("yp must be a real vector"));
1439 :     SET_SLOT(val, Matrix_ySym, duplicate(yp));
1440 :     if (!isMatrix(Xp) || !isReal(Xp))
1441 :     error(_("Xp must be a real matrix"));
1442 :     xdims = INTEGER(getAttrib(Xp, R_DimSymbol));
1443 :     if (xdims[0] != nobs) error(_("Xp must have %d rows"), nobs);
1444 :     p = xdims[1];
1445 :     xnms = VECTOR_ELT(getAttrib(Xp, R_DimNamesSymbol), 1);
1446 :     SET_SLOT(val, Matrix_XSym, duplicate(Xp));
1447 :     if (!isNewList(fl) || nf < 1) error(_("fl must be a nonempty list"));
1448 :     for (i = 0; i < nf; i++) {
1449 :     SEXP fli = VECTOR_ELT(fl, i);
1450 :     if (!isFactor(fli) || LENGTH(fli) != nobs)
1451 :     error(_("fl[[%d] must be a factor of length %d"), i+1, nobs);
1452 :     }
1453 :     SET_SLOT(val, Matrix_flistSym, duplicate(fl));
1454 :     if (!isString(method) || LENGTH(method) != 1)
1455 :     error(_("method must be a character vector of length 1"));
1456 :     SET_SLOT(val, Matrix_methodSym, duplicate(method));
1457 :     if (!isLogical(useS) || LENGTH(useS) != 1)
1458 :     error(_("useS must be a logical vector of length 1"));
1459 :     SET_SLOT(val, Matrix_useScaleSym, duplicate(useS));
1460 :     if (!isNewList(cnames) || LENGTH(cnames) != nf + 1)
1461 :     error(_("cnames must be a list of length %d"), nf + 1);
1462 :     SET_SLOT(val, Matrix_cnamesSym, duplicate(cnames));
1463 :     if (!isInteger(ncp) || LENGTH(ncp) != nf)
1464 :     error(_("ncp must be an integer vector of length %d"), nf);
1465 :     SET_SLOT(val, Matrix_callSym, duplicate(call));
1466 :     SET_SLOT(val, Matrix_familySym, duplicate(family));
1467 :     SET_SLOT(val, Matrix_ncSym, duplicate(ncp));
1468 :     Gp = INTEGER(ALLOC_SLOT(val, Matrix_GpSym, INTSXP, nf + 1));
1469 :     Gp[0] = 0;
1470 :     if (!isNewList(fl) || nf < 1) error(_("fl must be a nonempty list"));
1471 :     for (i = 0; i < nf; i++) {
1472 :     SEXP fli = VECTOR_ELT(fl, i);
1473 :     if (!isFactor(fli) || LENGTH(fli) != nobs)
1474 :     error(_("fl[[%d] must be a factor of length %d"), i+1, nobs);
1475 : bates 812
1476 : bates 1143 }
1477 :     SET_SLOT(val, Matrix_ZtSym, duplicate(ZZt));
1478 :     Zt = as_cholmod_sparse(GET_SLOT(val, Matrix_ZtSym));
1479 :     /* analyze Zt */
1480 :     q = Zt->nrow;
1481 :     i = c.supernodal;
1482 :     c.supernodal = CHOLMOD_SUPERNODAL; /* force a supernodal decomposition */
1483 :     F = cholmod_analyze(Zt, &c);
1484 :     c.supernodal = i;
1485 :     ts1 = cholmod_aat(Zt, (int*)NULL/* fset */,(size_t)0,
1486 :     CHOLMOD_PATTERN, &c);
1487 :     ts2 = cholmod_copy(ts1, 1/*upper triangle*/, CHOLMOD_PATTERN, &c);
1488 :     SET_SLOT(val, Matrix_ZtZSym,
1489 :     alloc_dsCMatrix(q, cholmod_nnz(ts2, &c), "U", R_NilValue,
1490 :     R_NilValue));
1491 :     cholmod_free_sparse(&ts1, &c); cholmod_free_sparse(&ts2, &c);
1492 :     /* allocate other slots */
1493 :     SET_SLOT(val, Matrix_devianceSym, Matrix_make_named(REALSXP, devnms));
1494 :     SET_SLOT(val, Matrix_devCompSym, Matrix_make_named(REALSXP, devCmpnms));
1495 :     SET_SLOT(val, Matrix_statusSym, Matrix_make_named(LGLSXP, statusnms));
1496 :     AZERO(LOGICAL(GET_SLOT(val, Matrix_statusSym)), 4);
1497 :     dcmp = REAL(GET_SLOT(val, Matrix_devCompSym));
1498 :     AZERO(dcmp, 7); /* cosmetic */
1499 :     dcmp[0] = (double) nobs;
1500 :     dcmp[1] = (double) p;
1501 :     /* allocate and populate list slots */
1502 :     Omega = ALLOC_SLOT(val, Matrix_OmegaSym, VECSXP, nf);
1503 :     bVar = ALLOC_SLOT(val, Matrix_bVarSym, VECSXP, nf);
1504 :     gradComp = ALLOC_SLOT(val, Matrix_gradCompSym, VECSXP, nf);
1505 :     setAttrib(Omega, R_NamesSymbol, duplicate(fnms));
1506 :     setAttrib(bVar, R_NamesSymbol, duplicate(fnms));
1507 :     setAttrib(gradComp, R_NamesSymbol, duplicate(fnms));
1508 :     for (i = 0; i < nf; i++) {
1509 :     int nci = nc[i];
1510 :     int nlev = LENGTH(getAttrib(VECTOR_ELT(fl, i), R_LevelsSymbol));
1511 :     SET_VECTOR_ELT(Omega, i,
1512 :     alloc_dpoMatrix(nci, "U",
1513 :     VECTOR_ELT(cnames, i),
1514 :     VECTOR_ELT(cnames, i)));
1515 :     SET_VECTOR_ELT(bVar, i, alloc3Darray(REALSXP, nci, nci, nlev));
1516 :     SET_VECTOR_ELT(gradComp, i, alloc3Darray(REALSXP, nci, nci, 4));
1517 :     Gp[i + 1] = Gp[i] + nlev * nci;
1518 :     }
1519 :     /* create ZtX, RZX, XtX, RXX */
1520 :     SET_SLOT(val, Matrix_ZtXSym, alloc_dgeMatrix(q, p, R_NilValue, xnms));
1521 :     SET_SLOT(val, Matrix_RZXSym, alloc_dgeMatrix(q, p, R_NilValue, xnms));
1522 :     SET_SLOT(val, Matrix_XtXSym, alloc_dpoMatrix(p, "U", xnms, xnms));
1523 :     SET_SLOT(val, Matrix_RXXSym, alloc_dtrMatrix(p, "U", "N", xnms, xnms));
1524 :     SET_SLOT(val, Matrix_ZtySym, allocVector(REALSXP, q));
1525 :     SET_SLOT(val, Matrix_rZySym, allocVector(REALSXP, q));
1526 :     SET_SLOT(val, Matrix_XtySym, allocVector(REALSXP, p));
1527 :     SET_SLOT(val, Matrix_rXySym, allocVector(REALSXP, p));
1528 : bates 1156 /* create weights and working residuals */
1529 :     wts = REAL(ALLOC_SLOT(val, Matrix_wtsSym, REALSXP, nobs));
1530 :     wrkres = REAL(ALLOC_SLOT(val, Matrix_wrkresSym, REALSXP, nobs));
1531 :     y = REAL(GET_SLOT(val, Matrix_ySym));
1532 :     for (i = 0; i < nobs; i++) {wts[i] = 1.; wrkres[i] = y[i];}
1533 : bates 1143 internal_mer_update_ZXy(val, (int*)(F->Perm));
1534 :     Free(Zt);
1535 :     /* secondary slots */
1536 :     SET_SLOT(val, Matrix_ranefSym, allocVector(REALSXP, q));
1537 :     SET_SLOT(val, Matrix_fixefSym, allocVector(REALSXP, p));
1538 :     SET_SLOT(val, Matrix_RZXinvSym, alloc_dgeMatrix(q, p, R_NilValue, xnms));
1539 :     /* initialize */
1540 :     mer_initial(val);
1541 :     /* The next calls are simply to set up the L slot. At present the
1542 :     * factor F is a symbolic factor. We need to force it to become
1543 :     * numeric before allocating the L slot in the object. */
1544 :     internal_mer_Zfactor(val, F);
1545 : bates 1156 /* One side-effect of this call is to set the status as
1546 : bates 1143 * factored. We need to reset it */
1547 :     LOGICAL(GET_SLOT(val, Matrix_statusSym))[0] = 0;
1548 :     /* Create the dCHMfactor object and store it in the L slot. This
1549 :     * also frees the storage. */
1550 :     SET_SLOT(val, Matrix_LSym, chm_factor_to_SEXP(F, 1));
1551 :     /* OK, done now. */
1552 : bates 812 UNPROTECT(1);
1553 : bates 1143 return val;
1554 : bates 812 }
1555 : bates 815
1556 : bates 1143 /**
1557 :     * Extract parameters from the Omega matrices. These aren't
1558 :     * "coefficients" but the extractor is called coef for historical
1559 :     * reasons. Within each group these values are in the order of the
1560 :     * diagonal entries first then the strict upper triangle in row
1561 :     * order.
1562 : bates 815 *
1563 : bates 1143 * The parameters can be returned in three forms:
1564 :     * 0 - nonlinearly constrained - elements of the relative precision matrix
1565 :     * 1 - unconstrained - from the LDL' decomposition - logarithms of
1566 :     * the diagonal elements of D
1567 :     * 2 - box constrained - also from the LDL' decomposition - inverses
1568 :     * of the diagonal elements of D
1569 :     *
1570 : bates 815 * @param x pointer to an mer object
1571 : bates 1143 * @param pType pointer to an integer scalar indicating the form of the
1572 :     * parameters to be returned.
1573 :     *
1574 :     * @return numeric vector of the values in the upper triangles of the
1575 :     * Omega matrices
1576 : bates 815 */
1577 : bates 1143 SEXP mer_coef(SEXP x, SEXP pType)
1578 : bates 815 {
1579 : bates 1143 int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1580 :     nf = LENGTH(GET_SLOT(x, Matrix_OmegaSym));
1581 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1582 : bates 815
1583 : bates 1143 internal_mer_coef(x, asInteger(pType), REAL(val));
1584 :     UNPROTECT(1);
1585 :     return val;
1586 :     }
1587 : bates 818
1588 : bates 1143 /**
1589 :     * Assign the upper triangles of the Omega matrices according to a
1590 :     * vector of parameters.
1591 :     *
1592 :     * @param x pointer to an lme object
1593 :     * @param coef pointer to an numeric vector of appropriate length
1594 :     * @param pType pointer to an integer scalar
1595 :     *
1596 :     * @return R_NilValue
1597 :     */
1598 :     SEXP mer_coefGets(SEXP x, SEXP coef, SEXP pType)
1599 :     {
1600 :     int clen = coef_length(LENGTH(GET_SLOT(x, Matrix_flistSym)),
1601 :     INTEGER(GET_SLOT(x, Matrix_ncSym)));
1602 :     if (LENGTH(coef) != clen || !isReal(coef))
1603 :     error(_("coef must be a numeric vector of length %d"), clen);
1604 :     internal_mer_coefGets(x, REAL(coef), asInteger(pType));
1605 :     return x;
1606 :     }
1607 : bates 818
1608 : bates 1143 /**
1609 :     * Return L as a dtCMatrix object
1610 :     *
1611 :     * @param x pointer to an mer object
1612 :     *
1613 :     * @return L as an dtCMatrix object
1614 :     */
1615 :     SEXP mer_dtCMatrix(SEXP x)
1616 :     {
1617 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym)), *Lcp;
1618 :     cholmod_sparse *Lm;
1619 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
1620 :     int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
1621 :     nz, q;
1622 : bates 818
1623 : bates 1143 dims[0] = dims[1] = q = (int)(L->n);
1624 :     Lcp = cholmod_copy_factor(L, &c); Free(L); /* next call changes Lcp */
1625 :     Lm = cholmod_factor_to_sparse(Lcp, &c); cholmod_free_factor(&Lcp, &c);
1626 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
1627 :     SET_SLOT(ans, Matrix_diagSym, mkString("N"));
1628 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, q + 1)),
1629 :     (int*) Lm->p, q + 1);
1630 :     nz = ((int*)(Lm->p))[q];
1631 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
1632 :     (int*) Lm->i, nz);
1633 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)),
1634 :     (double*) Lm->x, nz);
1635 :     cholmod_free_sparse(&Lm, &c);
1636 :     UNPROTECT(1);
1637 :     return ans;
1638 : bates 815 }
1639 :    
1640 :     /**
1641 : bates 1143 * Return L^{-1} as a dtCMatrix object
1642 : bates 815 *
1643 : bates 1143 * @param x pointer to an mer object
1644 : bates 815 *
1645 : bates 1143 * @return L^{-1} as an dtCMatrix object
1646 : bates 815 */
1647 : bates 1143 SEXP mer_dtCMatrix_inv(SEXP x)
1648 : bates 815 {
1649 : bates 1143 cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
1650 :     cholmod_sparse
1651 :     *b = cholmod_allocate_sparse(L->n, L->n, L->n, 1, 1,
1652 :     0, CHOLMOD_REAL, &c),
1653 :     *Linv;
1654 :     double *bx = (double *)(b->x);
1655 :     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtCMatrix")));
1656 :     int *bi = (int *) (b->i), *bp = (int *) (b->p),
1657 :     *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),
1658 :     j, nz, q;
1659 : bates 815
1660 : bates 1143 dims[0] = dims[1] = q = (int)(L->n);
1661 :     for (j = 0; j < q; j++) {
1662 :     bp[j] = bi[j] = j;
1663 :     bx[j] = 1;
1664 : bates 815 }
1665 : bates 1143 bp[q] = q;
1666 :     Linv = cholmod_spsolve(CHOLMOD_L, L, b, &c);
1667 :     cholmod_free_sparse(&b, &c);
1668 :     SET_SLOT(ans, Matrix_uploSym, mkString("L"));
1669 :     SET_SLOT(ans, Matrix_diagSym, mkString("N"));
1670 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, q + 1)),
1671 :     (int *) Linv->p, q + 1);
1672 :     nz = ((int *)(Linv->p))[q];
1673 :     Memcpy(INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz)),
1674 :     (int *) Linv->i, nz);
1675 :     Memcpy(REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz)),
1676 :     (double *) Linv->x, nz);
1677 :     cholmod_free_sparse(&Linv, &c);
1678 : bates 815 UNPROTECT(1);
1679 :     return ans;
1680 :     }
1681 : bates 879
1682 : bates 1143 /**
1683 :     * Create and factor Z'Z+Omega. Also create RZX and RXX, the deviance
1684 :     * components, and the value of the deviance for both ML and REML.
1685 :     *
1686 :     * @param x pointer to an lmer object
1687 :     *
1688 :     * @return NULL
1689 :     */
1690 :     SEXP mer_factor(SEXP x)
1691 :     {
1692 :     int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1693 :     if (!status[0]) {
1694 :     SEXP rXyP = GET_SLOT(x, Matrix_rXySym),
1695 :     rZyP = GET_SLOT(x, Matrix_rZySym);
1696 :     int i, info, ione = 1, p = LENGTH(rXyP), q = LENGTH(rZyP);
1697 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
1698 :     double *RXX = REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym), Matrix_xSym)),
1699 :     *RZX = REAL(GET_SLOT(GET_SLOT(x, Matrix_RZXSym), Matrix_xSym)),
1700 :     *rXy = REAL(rXyP), *rZy = REAL(rZyP),
1701 :     *dcmp = REAL(GET_SLOT(x, Matrix_devCompSym)),
1702 :     *dev = REAL(GET_SLOT(x, Matrix_devianceSym)),
1703 :     one[2] = {1, 0}, m1[2] = {-1, 0};
1704 :     double nml = dcmp[0], nreml = dcmp[0] - dcmp[1];
1705 :    
1706 :     /* Inflate Z'Z to Z'Z+Omega and factor to form L. Form RZX and
1707 :     * rZy. Update status flags, dcmp[4] and dcmp[5]. */
1708 :     internal_mer_Zfactor(x, L);
1709 :     /* downdate XtX and factor */
1710 :     Memcpy(RXX, REAL(GET_SLOT(GET_SLOT(x, Matrix_XtXSym), Matrix_xSym)), p * p);
1711 :     F77_CALL(dsyrk)("U", "T", &p, &q, m1, RZX, &q, one, RXX, &p);
1712 :     F77_CALL(dpotrf)("U", &p, RXX, &p, &info);
1713 :     if (info) {
1714 :     error(_("Leading minor of order %d in downdated X'X is not positive definite"),
1715 :     info);
1716 :     dcmp[3] = dcmp[6] = dev[0] = dev[1] = NA_REAL;
1717 :     } else {
1718 :     for (dcmp[6] = 0, i = 0; i < p; i++) /* 2 * logDet(RXX) */
1719 :     dcmp[6] += 2. * log(RXX[i * (p + 1)]);
1720 :     /* solve for rXy and ryy^2 */
1721 :     Memcpy(rXy, REAL(GET_SLOT(x, Matrix_XtySym)), p);
1722 :     F77_CALL(dgemv)("T", &q, &p, m1, RZX, &q, rZy, &ione, one, rXy, &ione);
1723 :     F77_CALL(dtrsv)("U", "T", "N", &p, RXX, &p, rXy, &ione);
1724 :     dcmp[3] = log(dcmp[2] /* dcmp[3] = log(ryy^2); dcmp[2] = y'y; */
1725 :     - F77_CALL(ddot)(&p, rXy, &ione, rXy, &ione)
1726 :     - F77_CALL(ddot)(&q, rZy, &ione, rZy, &ione));
1727 :     /* evaluate ML and REML deviance */
1728 :     dev[0] = dcmp[4] - dcmp[5] +
1729 :     nml*(1.+dcmp[3]+log(2.*PI/nml));
1730 :     dev[1] = dcmp[4] - dcmp[5] + dcmp[6] +
1731 :     nreml*(1.+dcmp[3]+log(2.*PI/nreml));
1732 :     }
1733 :     Free(L);
1734 :     }
1735 :     return R_NilValue;
1736 :     }
1737 :    
1738 : bates 879 /**
1739 : bates 1143 * Return the fitted values as an SEXP
1740 : bates 879 *
1741 : bates 1143 * @param x pointer to an mer object
1742 :     * @param useFe pointer to a logical scalar indicating if the fixed
1743 :     * effects should be used
1744 :     * @param useRe pointer to a logical scalar indicating if the random
1745 :     * effects should be used
1746 : bates 879 *
1747 : bates 1143 * @return pointer to a numeric array of fitted values
1748 : bates 879 */
1749 : bates 1143
1750 : bates 1162 SEXP mer_fitted(SEXP x)
1751 : bates 879 {
1752 : bates 1143 int n = LENGTH(GET_SLOT(x, Matrix_ySym));
1753 :     SEXP ans = PROTECT(allocVector(REALSXP, n));
1754 : bates 879
1755 : bates 1162 internal_mer_fitted(x, (double*) NULL, REAL(ans));
1756 : bates 879 UNPROTECT(1);
1757 :     return ans;
1758 :     }
1759 : bates 888
1760 : bates 1143 /**
1761 :     * Extract the conditional estimates of the fixed effects
1762 :     *
1763 : bates 888 * @param x Pointer to an mer object
1764 : bates 1143 *
1765 :     * @return a numeric vector containing the conditional estimates of
1766 :     * the fixed effects
1767 : bates 888 */
1768 : bates 1143 SEXP mer_fixef(SEXP x)
1769 : bates 888 {
1770 : bates 1143 int nf = LENGTH(GET_SLOT(x, Matrix_OmegaSym));
1771 :     SEXP ans;
1772 : bates 888
1773 : bates 1143 mer_secondary(x);
1774 :     ans = PROTECT(duplicate(GET_SLOT(x, Matrix_fixefSym)));
1775 :     setAttrib(ans, R_NamesSymbol,
1776 :     duplicate(VECTOR_ELT(GET_SLOT(x, Matrix_cnamesSym), nf)));
1777 :     UNPROTECT(1);
1778 :     return ans;
1779 : bates 888 }
1780 : bates 940
1781 : bates 1143 /**
1782 :     * Fill in the gradComp and bVar slots. Each component in the gradComp slot
1783 :     * consists of four symmetric matrices used to generate the gradient or the ECME
1784 :     * step. They are
1785 :     * 1) -m_i\bOmega_i^{-1}
1786 :     * 2) \bB_i\bB_i\trans
1787 :     * 3) \tr\left[\der_{\bOmega_i}\bOmega\left(\bZ\trans\bZ+\bOmega\right)\inv\right]
1788 :     * 4) The term added to 3) to get \tr\left[\der_{\bOmega_i}\bOmega\vb\right]
1789 :     *
1790 :     * @param x pointer to an mer object
1791 :     * @param val pointer to a list of matrices of the correct sizes
1792 :     *
1793 :     * @return NULL
1794 :     */
1795 :     SEXP mer_gradComp(SEXP x)
1796 : bates 940 {
1797 : bates 1143 int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
1798 : bates 940
1799 : bates 1143 mer_secondary(x);
1800 :     if (!status[2]) {
1801 :     SEXP bVarP = GET_SLOT(x, Matrix_bVarSym),
1802 :     OmegaP = GET_SLOT(x, Matrix_OmegaSym),
1803 :     RZXP = GET_SLOT(x, Matrix_RZXSym),
1804 :     gradComp = GET_SLOT(x, Matrix_gradCompSym),
1805 :     ranefP = GET_SLOT(x, Matrix_ranefSym);
1806 :     int q = LENGTH(ranefP), p = LENGTH(GET_SLOT(x, Matrix_rXySym));
1807 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
1808 :     cholmod_dense *RZX = as_cholmod_dense(RZXP), *tmp1;
1809 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1810 :     *Perm = (int *)(L->Perm),
1811 :     *iperm = Calloc(q, int),
1812 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1813 :     i, j, k, nf = length(OmegaP);
1814 :     double *b = REAL(GET_SLOT(x, Matrix_ranefSym)),
1815 :     *RZXinv = REAL(GET_SLOT(GET_SLOT(x, Matrix_RZXinvSym), Matrix_xSym)),
1816 :     m1[] = {-1, 0}, alpha;
1817 :    
1818 :     alpha = 1./internal_mer_sigma(x, -1);
1819 :     alpha = alpha * alpha;
1820 :    
1821 :     for (j = 0; j < q; j++) iperm[Perm[j]] = j; /* create the inverse permutation */
1822 :     /* form RZXinv, the section of R^{-1} from RZX */
1823 :     tmp1 = cholmod_solve(CHOLMOD_Lt, L, RZX, &c); Free(RZX);
1824 :     /* copy columns of tmp1 to RZXinv applying the inverse permutation */
1825 :     for (j = 0; j < p; j++) {
1826 :     double *dest = RZXinv + j * q, *src = ((double*)(tmp1->x)) + j * q;
1827 :     for (i = 0; i < q; i++) dest[i] = src[iperm[i]];
1828 :     }
1829 :     cholmod_free_dense(&tmp1, &c);
1830 :     F77_CALL(dtrsm)("R", "U", "N", "N", &q, &p, m1,
1831 :     REAL(GET_SLOT(GET_SLOT(x, Matrix_RXXSym), Matrix_xSym)),
1832 :     &p, RZXinv, &q);
1833 : bates 940
1834 : bates 1143 internal_mer_bVar(x);
1835 :     for (i = 0; i < nf; i++) {
1836 :     int nci = nc[i], RZXrows = Gp[i + 1] - Gp[i];
1837 :     int ncisq = nci * nci, nlev = RZXrows/nci;
1838 :     double *bVi = REAL(VECTOR_ELT(bVarP, i)),
1839 :     *bi = b + Gp[i], *mm = REAL(VECTOR_ELT(gradComp, i)),
1840 :     *tmp = Memcpy(Calloc(ncisq, double),
1841 :     REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(OmegaP, i)),
1842 :     Matrix_xSym)), ncisq),
1843 :     *RZXi = RZXinv + Gp[i],
1844 :     dlev = (double) nlev,
1845 :     one[] = {1,0}, zero[] = {0,0};
1846 : bates 940
1847 : bates 1143 if (nci == 1) {
1848 :     int ione = 1;
1849 :     mm[0] = ((double) nlev)/(tmp[0] * tmp[0]);
1850 :     mm[1] = alpha * F77_CALL(ddot)(&nlev, bi, &ione, bi, &ione);
1851 :     mm[2] = 0.;
1852 :     for (k = 0; k < nlev; k++) mm[2] += bVi[k];
1853 :     mm[3] = 0.;
1854 :     for (j = 0; j < p; j++) {
1855 :     mm[3] += F77_CALL(ddot)(&RZXrows, RZXi + j * q, &ione,
1856 :     RZXi + j * q, &ione);
1857 :     }
1858 :     } else {
1859 :     AZERO(mm, 4 * ncisq);
1860 :     F77_CALL(dtrtri)("U", "N", &nci, tmp, &nci, &j);
1861 :     if (j)
1862 :     error(_("Omega[[%d]] is not positive definite"), i + 1);
1863 :     F77_CALL(dsyrk)("U", "N", &nci, &nci, &dlev, tmp, &nci,
1864 :     zero, mm, &nci);
1865 :     mm += ncisq; /* \bB_i term */
1866 :     F77_CALL(dsyrk)("U", "N", &nci, &nlev, &alpha, bi, &nci,
1867 :     zero, mm, &nci);
1868 :     mm += ncisq; /* Sum of diagonal blocks of the inverse
1869 :     * (Z'Z+Omega)^{-1} */
1870 :     for (j = 0; j < ncisq; j++) {
1871 :     for (k = 0; k < nlev; k++) mm[j] += bVi[j + k*ncisq];
1872 :     }
1873 :     mm += ncisq; /* Extra term for \vb */
1874 :     for (j = 0; j < p; j++) {
1875 :     F77_CALL(dsyrk)("U", "N", &nci, &nlev, one,
1876 :     RZXi + j * q, &nci,
1877 :     one, mm, &nci);
1878 :     }
1879 : bates 963 }
1880 : bates 1143 Free(tmp);
1881 : bates 940 }
1882 : bates 1143 Free(iperm); Free(L);
1883 :     status[2] = 1; status[3] = 0;
1884 : bates 940 }
1885 : bates 1143 return R_NilValue;
1886 : bates 940 }
1887 :    
1888 : bates 1143 /**
1889 :     * Evaluate the gradient vector
1890 :     *
1891 :     * @param x Pointer to an lmer object
1892 :     * @param pType Pointer to an integer indicator of the
1893 :     * parameterization being used
1894 :     *
1895 :     * @return pointer to a gradient vector
1896 :     */
1897 :     SEXP mer_gradient(SEXP x, SEXP pType)
1898 : bates 967 {
1899 : bates 1143 SEXP Omega = GET_SLOT(x, Matrix_OmegaSym);
1900 :     SEXP gradComp = GET_SLOT(x, Matrix_gradCompSym);
1901 :     int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1902 :     dind, i, ifour = 4, ione = 1, nf = length(Omega),
1903 :     odind, ptyp = asInteger(pType);
1904 :     int REML =
1905 :     !strcmp(CHAR(asChar(GET_SLOT(x, Matrix_methodSym))), "REML");
1906 :     SEXP val = PROTECT(allocVector(REALSXP, coef_length(nf, nc)));
1907 :     double cc[] = {-1, 1, 1, REML ? 1 : 0},
1908 :     one = 1.0, zero = 0.0;
1909 : bates 967
1910 : bates 1143 mer_gradComp(x);
1911 :     dind = 0; /* index into val for diagonals */
1912 :     for (i = 0; i < nf; i++) {
1913 :     SEXP Omgi = VECTOR_ELT(Omega, i);
1914 :     int nci = nc[i], ncisqr = nci * nci;
1915 :     double *tmp = Calloc(ncisqr, double);
1916 : bates 967
1917 : bates 1143 F77_CALL(dgemm)("N", "N", &ncisqr, &ione, &ifour, &one,
1918 :     REAL(VECTOR_ELT(gradComp, i)), &ncisqr,
1919 :     cc, &ifour, &zero, tmp, &ncisqr);
1920 :     if (nci == 1) {
1921 :     double omg = *REAL(GET_SLOT(Omgi, Matrix_xSym));
1922 :     REAL(val)[dind++] =
1923 :     (ptyp?((ptyp == 1)? omg : -omg * omg) : 1) * tmp[0];
1924 :     } else {
1925 :     int ii, j, ncip1 = nci + 1;
1926 : bates 967
1927 : bates 1143 odind = dind + nci; /* index into val for off-diagonals */
1928 :     if (ptyp) {
1929 :     double *chol = REAL(GET_SLOT(dpoMatrix_chol(Omgi), Matrix_xSym)),
1930 :     *tmp2 = Calloc(ncisqr, double);
1931 : bates 967
1932 : bates 1143 /* Overwrite the gradient with respect to positions in
1933 :     * Omega[[i]] by the gradient with respect to the
1934 :     * unconstrained parameters.*/
1935 : bates 967
1936 : bates 1143 /* tmp2 := chol %*% tmp using only upper triangle of tmp */
1937 :     F77_CALL(dsymm)("R", "U", &nci, &nci, &one, tmp, &nci,
1938 :     chol, &nci, &zero, tmp2, &nci);
1939 :     /* full symmetric product gives diagonals */
1940 :     F77_CALL(dtrmm)("R", "U", "T", "N", &nci, &nci, &one, chol,
1941 :     &nci, Memcpy(tmp, tmp2, ncisqr), &nci);
1942 :     /* overwrite upper triangle with gradients for L' */
1943 :     for (ii = 1; ii < nci; ii++) {
1944 :     for (j = 0; j < ii; j++) {
1945 :     tmp[j + ii*nci] = chol[j*ncip1] * tmp2[j + ii*nci];
1946 :     tmp[ii + j*nci] = 0.;
1947 :     }
1948 :     }
1949 :     if (ptyp > 1)
1950 :     for (ii = 0; ii < nci; ii++) {
1951 :     int ind = ii * ncip1;
1952 :     double sqrtd = chol[ind];
1953 :     tmp[ind] *= -(sqrtd*sqrtd);
1954 :     }
1955 :     }
1956 :     for (j = 0; j < nci; j++) {
1957 :     REAL(val)[dind + j] = tmp[j * ncip1];
1958 :     for (ii = 0; ii < j; ii++) /* offdiagonals count twice */
1959 :     REAL(val)[odind++] = 2. * tmp[ii + j * nci];
1960 :     }
1961 :     dind = odind;
1962 :     }
1963 :     Free(tmp);
1964 :     }
1965 : bates 967 UNPROTECT(1);
1966 : bates 1143 return val;
1967 : bates 967 }
1968 :    
1969 : bates 940 /**
1970 : bates 978 * Create and insert initial values for Omega.
1971 :     *
1972 : bates 1143 * @param x pointer to an mer object
1973 : bates 978 *
1974 :     * @return NULL
1975 :     */
1976 : bates 1143 SEXP mer_initial(SEXP x)
1977 : bates 978 {
1978 :     SEXP Omg = GET_SLOT(x, Matrix_OmegaSym),
1979 :     ZtZ = GET_SLOT(x, Matrix_ZtZSym);
1980 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
1981 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
1982 :     *p = INTEGER(GET_SLOT(ZtZ, Matrix_pSym)),
1983 :     i, nf = length(Omg);
1984 :     double *xx = REAL(GET_SLOT(ZtZ, Matrix_xSym));
1985 :    
1986 :     for (i = 0; i < nf; i++) {
1987 : bates 1143 SEXP Omgi = VECTOR_ELT(Omg, i);
1988 :     double *omgi = REAL(GET_SLOT(Omgi, Matrix_xSym));
1989 : bates 978 int bb = Gp[i], j, k, nci = nc[i];
1990 :     int ncip1 = nci + 1, nlev = (Gp[i + 1] - bb)/nci;
1991 :    
1992 :     AZERO(omgi, nci * nci);
1993 :     for (j = 0; j < nlev; j++) {
1994 :     int base = bb + j * nci;
1995 :     for (k = 0; k < nci; k++)
1996 : bates 1143 /* add the last element in the column */
1997 : bates 978 omgi[k * ncip1] += xx[p[base + k + 1] - 1];
1998 :     }
1999 :     for (k = 0; k < nci; k++) omgi[k * ncip1] *= 0.375/nlev;
2000 : bates 1143 SET_SLOT(Omgi, Matrix_factorSym, allocVector(VECSXP, 0));
2001 :     dpoMatrix_chol(Omgi);
2002 : bates 978 }
2003 : bates 1143 flag_not_factored(x);
2004 : bates 978 return R_NilValue;
2005 :     }
2006 :    
2007 :     /**
2008 : bates 1143 * Extract the conditional modes of the random effects.
2009 : bates 940 *
2010 : bates 1143 * @param x Pointer to an mer object
2011 : bates 940 *
2012 : bates 1143 * @return a list of matrices containing the conditional modes of the
2013 :     * random effects
2014 : bates 940 */
2015 : bates 1143 SEXP mer_ranef(SEXP x)
2016 : bates 940 {
2017 : bates 1143 SEXP cnames = GET_SLOT(x, Matrix_cnamesSym),
2018 :     flist = GET_SLOT(x, Matrix_flistSym);
2019 :     int *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
2020 :     *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
2021 :     i, ii, jj,
2022 :     nf = length(flist);
2023 :     SEXP val = PROTECT(allocVector(VECSXP, nf));
2024 :     double *b = REAL(GET_SLOT(x, Matrix_ranefSym));
2025 : bates 963
2026 : bates 1143 mer_secondary(x);
2027 :     setAttrib(val, R_NamesSymbol,
2028 :     duplicate(getAttrib(flist, R_NamesSymbol)));
2029 :     for (i = 0; i < nf; i++) {
2030 :     SEXP nms, rnms = getAttrib(VECTOR_ELT(flist, i), R_LevelsSymbol);
2031 :     int nci = nc[i], mi = length(rnms);
2032 :     double *bi = b + Gp[i], *mm;
2033 :    
2034 :     SET_VECTOR_ELT(val, i, allocMatrix(REALSXP, mi, nci));
2035 :     setAttrib(VECTOR_ELT(val, i), R_DimNamesSymbol, allocVector(VECSXP, 2));
2036 :     nms = getAttrib(VECTOR_ELT(val, i), R_DimNamesSymbol);
2037 :     SET_VECTOR_ELT(nms, 0, duplicate(rnms));
2038 :     SET_VECTOR_ELT(nms, 1, duplicate(VECTOR_ELT(cnames, i)));
2039 :     mm = REAL(VECTOR_ELT(val, i));
2040 :     for (jj = 0; jj < nci; jj++)
2041 :     for(ii = 0; ii < mi; ii++)
2042 :     mm[ii + jj * mi] = bi[jj + ii * nci];
2043 :     }
2044 : bates 940 UNPROTECT(1);
2045 :     return val;
2046 :     }
2047 : bates 967
2048 : bates 1143 /**
2049 :     * Update the secondary slots - fixef and ranef
2050 :     *
2051 :     * @param x pointer to a mer object
2052 :     *
2053 :     */
2054 :     SEXP mer_secondary(SEXP x)
2055 : bates 977 {
2056 : bates 1143 int *status = LOGICAL(GET_SLOT(x, Matrix_statusSym));
2057 : bates 977
2058 : bates 1143 mer_factor(x);
2059 :     if (!status[1]) {
2060 :     internal_mer_fixef(x);
2061 :     internal_mer_ranef(x);
2062 : bates 978 }
2063 : bates 1143 return R_NilValue;
2064 : bates 978 }
2065 :    
2066 : bates 967 /**
2067 : bates 1143 * Extract the ML or REML conditional estimate of sigma
2068 : bates 967 *
2069 : bates 1143 * @param x pointer to an mer object
2070 :     * @param REML logical scalar - TRUE if REML estimates are requested
2071 : bates 967 *
2072 : bates 1143 * @return pointer to a numeric scalar
2073 : bates 967 */
2074 : bates 1143 SEXP mer_sigma(SEXP x, SEXP REML)
2075 : bates 967 {
2076 : bates 1143 return ScalarReal(
2077 :     internal_mer_sigma(x,
2078 :     (REML == R_NilValue) ? -1 :
2079 :     (asLogical(REML))));
2080 : bates 967 }
2081 : bates 978
2082 :     /**
2083 : bates 1143 * Simulate a set of linear predictors from the random effects part of
2084 :     * an mer object
2085 : bates 978 *
2086 : bates 1143 * @param x Pointer to an mer object
2087 :     * @param np Pointer to an integer giving the number of values to simulate
2088 :     * @param useScale Logical indicator of whether to use the scale
2089 : bates 978 *
2090 : bates 1143 * @return a matrix of simulated linear predictors
2091 : bates 978 */
2092 : bates 1143 SEXP mer_simulate(SEXP x, SEXP nsimP)
2093 : bates 978 {
2094 : bates 1143 int *nc = INTEGER(GET_SLOT(x, Matrix_ncSym)),
2095 :     *Gp = INTEGER(GET_SLOT(x, Matrix_GpSym)),
2096 :     REML = !strcmp(CHAR(asChar(GET_SLOT(x, Matrix_methodSym))),"REML"),
2097 :     i, ii, j, nsim = asInteger(nsimP),
2098 :     nf = LENGTH(GET_SLOT(x, Matrix_OmegaSym)),
2099 :     n = LENGTH(GET_SLOT(x, Matrix_ySym)),
2100 :     q = LENGTH(GET_SLOT(x, Matrix_ZtySym));
2101 :     SEXP ans = PROTECT(allocMatrix(REALSXP, n, nsim)),
2102 :     Omega = GET_SLOT(x, Matrix_OmegaSym);
2103 :     cholmod_dense *cha = as_cholmod_dense(ans),
2104 :     *chb = cholmod_allocate_dense(q, nsim, q, CHOLMOD_REAL, &c);
2105 :     double one[] = {1,0}, zero[] = {0,0},
2106 :     scale = (asLogical(GET_SLOT(x, Matrix_useScaleSym)) ?
2107 :     internal_mer_sigma(x, REML) : 1);
2108 :     cholmod_sparse *Zt = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtSym));
2109 :    
2110 :     GetRNGstate();
2111 :     for (ii = 0; ii < nsim; ii++) {
2112 :     for (i = 0; i < nf; i++) {
2113 :     int nci = nc[i], relen = Gp[i + 1] - Gp[i];
2114 :     int nlev = relen/nci;
2115 :     double *bi = (double *)(chb->x) + ii * q + Gp[i],
2116 :     *Rmat = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)),
2117 :     Matrix_xSym));
2118 : bates 978
2119 : bates 1143 for (j = 0; j < relen; j++) bi[j] = norm_rand();
2120 :     F77_CALL(dtrsm)("L", "U", "N", "N", &nci, &nlev, &scale,
2121 :     Rmat, &nci, bi, &nci);
2122 :     }
2123 :     }
2124 :     PutRNGstate();
2125 :    
2126 :     if (!cholmod_sdmult(Zt, 1, one, zero, chb, cha, &c))
2127 :     error(_("cholmod_sdmult failed"));
2128 :     cholmod_free_dense(&chb, &c);
2129 :     Free(Zt); Free(cha);
2130 : bates 978 UNPROTECT(1);
2131 :     return ans;
2132 :     }
2133 :    
2134 :     /**
2135 : bates 1143 * Externally callable version of internal_mer_update_ZXy
2136 : bates 978 *
2137 : bates 1143 * @param x pointer to an mer object
2138 : bates 978 *
2139 : bates 1143 * @return NULL
2140 : bates 978 */
2141 : bates 1143 SEXP mer_update_ZXy(SEXP x)
2142 : bates 978 {
2143 : bates 1143 internal_mer_update_ZXy(x,
2144 :     INTEGER(GET_SLOT(GET_SLOT(x, Matrix_LSym),
2145 :     Matrix_permSym)));
2146 :     return R_NilValue;
2147 : bates 978 }
2148 :    
2149 : bates 1143 /**
2150 :     * Update the y slot (and slots derived from it) in an mer object
2151 :     *
2152 :     * @param x pointer to an mer object
2153 :     * @param ynew pointer to a numeric vector of length n
2154 :     *
2155 :     * @return NULL
2156 :     */
2157 :     SEXP mer_update_y(SEXP x, SEXP ynew)
2158 : bates 978 {
2159 : bates 1143 SEXP y = GET_SLOT(x, Matrix_ySym),
2160 :     Xty = GET_SLOT(x, Matrix_XtySym),
2161 :     Zty = GET_SLOT(x, Matrix_ZtySym);
2162 :     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));
2163 :     int *perm = (int*)(L->Perm), i, ione = 1,
2164 :     n = LENGTH(y), p = LENGTH(Xty), q = LENGTH(Zty);
2165 :     cholmod_sparse *Zt = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtSym));
2166 :     cholmod_dense *td1, *yd = as_cholmod_dense(GET_SLOT(x, Matrix_ySym));
2167 :     double one[] = {1,0}, zero[] = {0,0};
2168 : bates 978
2169 : bates 1143 if (!isReal(ynew) || LENGTH(ynew) != n)
2170 :     error(_("ynew must be a numeric vector of length %d"), n);
2171 :     Memcpy(REAL(y), REAL(ynew), n);
2172 :     /* y'y */
2173 :     REAL(GET_SLOT(x, Matrix_devCompSym))[2] =
2174 :     F77_CALL(ddot)(&n, REAL(y), &ione, REAL(y), &ione);
2175 :     /* PZ'y into Zty */
2176 :     td1 = cholmod_allocate_dense(q, 1, q, CHOLMOD_REAL, &c);
2177 :     if (!cholmod_sdmult(Zt, 0, one, zero, yd, td1, &c))
2178 :     error(_("cholmod_sdmult failed"));
2179 :     for (i = 0; i < q; i++) REAL(Zty)[i] = ((double *)(td1->x))[perm[i]];
2180 :     cholmod_free_dense(&td1, &c); Free(yd); Free(Zt);
2181 :     /* Xty */
2182 :     F77_CALL(dgemv)("T", &n, &p, one, REAL(GET_SLOT(x, Matrix_XSym)),
2183 :     &n, REAL(y), &ione, zero, REAL(Xty), &ione);
2184 :     flag_not_factored(x);
2185 :     Free(L);
2186 :     return R_NilValue;
2187 : bates 978
2188 :     }
2189 :    
2190 :     /**
2191 : bates 1143 * Check validity of an mer object.
2192 : bates 978 *
2193 : bates 1143 * @param x Pointer to an mer object
2194 : bates 978 *
2195 : bates 1143 * @return TRUE if the object is a valid lmer object, else a string
2196 :     * describing the nature of the violation.
2197 : bates 978 */
2198 : bates 1143 SEXP mer_validate(SEXP x)
2199 : bates 978 {
2200 : bates 1143 SEXP
2201 :     /* ZZxP = GET_SLOT(x, Matrix_ZZxSym), */
2202 :     ZtXP = GET_SLOT(x, Matrix_ZtXSym),
2203 :     XtXP = GET_SLOT(x, Matrix_XtXSym),
2204 :     RZXP = GET_SLOT(x, Matrix_RZXSym),
2205 :     RXXP = GET_SLOT(x, Matrix_RXXSym)
2206 :     /* , cnames = GET_SLOT(x, Matrix_cnamesSym) */
2207 :     ;
2208 :     int *ZtXd = INTEGER(getAttrib(ZtXP, Matrix_DimSym)),
2209 :     *XtXd = INTEGER(getAttrib(XtXP, Matrix_DimSym));
2210 : bates 978
2211 : bates 1143 if (!match_mat_dims(ZtXd, INTEGER(getAttrib(RZXP, Matrix_DimSym))))
2212 :     return mkString(_("Dimensions of slots ZtX and RZX must match"));
2213 :     if (!match_mat_dims(XtXd, INTEGER(getAttrib(RXXP, Matrix_DimSym))))
2214 :     return mkString(_("Dimensions of slots XtX and RXX must match"));
2215 :     if (ZtXd[1] != XtXd[0] || XtXd[0] != XtXd[1])
2216 :     return mkString(_("Slot XtX must be a square matrix with same ncol as ZtX"));
2217 :     return ScalarLogical(1);
2218 : bates 978 }
2219 :    
2220 : bates 1143 /**
2221 :     * Create the sparse Zt matrix from a factor list and list of model matrices
2222 :     *
2223 :     * @param fl list of factors
2224 :     * @param Ztl list of transposes of model matrices
2225 :     *
2226 :     * @return a freshly created sparse Zt object
2227 :     */
2228 :     SEXP Zt_create(SEXP fl, SEXP Ztl)
2229 : bates 978 {
2230 : bates 1143 SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgCMatrix"))), fi, tmmat;
2231 :     int *dims, *p, *ii, i, nrtot = 0, nf = LENGTH(fl), nobs;
2232 :     int *Gp = Calloc(nf + 1, int), *nr = Calloc(nf, int),
2233 :     *offset = Calloc(nf + 1, int);
2234 :     double *x;
2235 :    
2236 :     if (!isNewList(fl) || nf < 1)
2237 :     error(_("fl must be a non-null list"));
2238 :     if (!isNewList(Ztl) || LENGTH(Ztl) != nf)
2239 :     error(_("Ztl must be a list of length %d"), nf);
2240 :     fi = VECTOR_ELT(fl, 0);
2241 :     nobs = LENGTH(fi);
2242 :     if (!isFactor(fi) || nobs < 1)
2243 :     error(_("fl[[1]] must be a non-empty factor"));
2244 :     offset[0] = Gp[0] = 0;
2245 :     for (i = 0; i < nf; i++) { /* check consistency and establish dimensions */
2246 :     fi = VECTOR_ELT(fl, i); /* grouping factor */
2247 :     if (!isFactor(fi) || LENGTH(fi) != nobs)
2248 :     error(_("fl[[%d]] must be a factor of length %d"), i + 1, nobs);
2249 :     tmmat = VECTOR_ELT(Ztl, i); /* transpose of model matrix */
2250 :     if (!isMatrix(tmmat) || !isReal(tmmat))
2251 :     error(_("Ztl[[%d]] must be real matrix"), i + 1);
2252 :     dims = INTEGER(getAttrib(tmmat, R_DimSymbol));
2253 :     if (dims[1] != nobs)
2254 :     error(_("Ztl[[%d]] must have %d columns"), i + 1, nobs);
2255 :     nrtot += (nr[i] = dims[0]);
2256 :     offset[i + 1] = offset[i] + nr[i];
2257 :     Gp[i + 1] = Gp[i] + dims[0] * LENGTH(getAttrib(fi, R_LevelsSymbol));
2258 :     }
2259 :     dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));
2260 :     dims[0] = Gp[nf]; dims[1] = nobs;
2261 :     p = INTEGER(ALLOC_SLOT(ans, Matrix_pSym, INTSXP, nobs + 1));
2262 :     ii = INTEGER(ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nrtot * nobs));
2263 :     x = REAL(ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nrtot * nobs));
2264 :     p[0] = 0; for(i = 0; i < nobs; i++) p[i + 1] = p[i] + nrtot;
2265 : bates 978
2266 : bates 1143 for (i = 0; i < nf; i++) { /* fill ans */
2267 :     int *vals = INTEGER(VECTOR_ELT(fl, i)), j;
2268 :     double *xvals = REAL(VECTOR_ELT(Ztl, i));
2269 : bates 978
2270 : bates 1143 for (j = 0; j < nobs; j++) {
2271 :     int jbase = Gp[i] + nr[i] * (vals[j] - 1), k;
2272 :     for (k = 0; k < nr[i]; k++) {
2273 :     int ind = j * nrtot + offset[i] + k;
2274 :     ii[ind] = jbase + k;
2275 :     x[ind] = xvals[j * nr[i] + k];
2276 : bates 978 }
2277 :     }
2278 :     }
2279 :    
2280 : bates 1143 Free(offset); Free(Gp); Free(nr);
2281 :     UNPROTECT(1);
2282 :     return ans;
2283 : bates 978 }
2284 : bates 1162
2285 :     #if 0
2286 :    
2287 :     /* Gauss-Hermite Quadrature x positions and weights */
2288 :     static const double
2289 :     GHQ_x1[1] = {0},
2290 :     GHQ_w1[1] = {1},
2291 :     GHQ_x2[1] = {1},
2292 :     GHQ_w2[1] = {0.5},
2293 :     GHQ_x3[2] = {1.7320507779261, 0},
2294 :     GHQ_w3[2] = {0.166666666666667, 0.666666666666667},
2295 :     GHQ_x4[2] = {2.3344141783872, 0.74196377160456},
2296 :     GHQ_w4[2] = {0.0458758533899086, 0.454124131589555},
2297 :     GHQ_x5[3] = {2.85696996497785, 1.35562615677371, 0},
2298 :     GHQ_w5[3] = {0.0112574109895360, 0.222075915334214,
2299 :     0.533333317311434},
2300 :     GHQ_x6[3] = {3.32425737665988, 1.88917584542184,
2301 :     0.61670657963811},
2302 :     GHQ_w6[3] = {0.00255578432527774, 0.0886157433798025,
2303 :     0.408828457274383},
2304 :     GHQ_x7[4] = {3.7504396535397, 2.36675937022918,
2305 :     1.15440537498316, 0},
2306 :     GHQ_w7[4] = {0.000548268839501628, 0.0307571230436095,
2307 :     0.240123171391455, 0.457142843409801},
2308 :     GHQ_x8[4] = {4.14454711519499, 2.80248581332504,
2309 :     1.63651901442728, 0.539079802125417},
2310 :     GHQ_w8[4] = {0.000112614534992306, 0.0096352198313359,
2311 :     0.117239904139746, 0.373012246473389},
2312 :     GHQ_x9[5] = {4.51274578616743, 3.20542894799789,
2313 :     2.07684794313409, 1.02325564627686, 0},
2314 :     GHQ_w9[5] = {2.23458433364535e-05, 0.00278914123744297,
2315 :     0.0499164052656755, 0.244097495561989,
2316 :     0.406349194142045},
2317 :     GHQ_x10[5] = {4.85946274516615, 3.58182342225163,
2318 :     2.48432579912153, 1.46598906930182,
2319 :     0.484935699216176},
2320 :     GHQ_w10[5] = {4.31065250122166e-06, 0.000758070911538954,
2321 :     0.0191115799266379, 0.135483698910192,
2322 :     0.344642324578594},
2323 :     GHQ_x11[6] = {5.18800113558601, 3.93616653976536,
2324 :     2.86512311160915, 1.87603498804787,
2325 :     0.928868981484148, 0},
2326 :     GHQ_w11[6] = {8.12184954622583e-07, 0.000195671924393029,
2327 :     0.0067202850336527, 0.066138744084179,
2328 :     0.242240292596812, 0.36940835831095};
2329 :    
2330 :     static const double
2331 :     *GHQ_x[12] = {(double *) NULL, GHQ_x1, GHQ_x2, GHQ_x3, GHQ_x4,
2332 :     GHQ_x5, GHQ_x6, GHQ_x7, GHQ_x8, GHQ_x9, GHQ_x10,
2333 :     GHQ_x11},
2334 :     *GHQ_w[12] = {(double *) NULL, GHQ_w1, GHQ_w2, GHQ_w3, GHQ_w4,
2335 :     GHQ_w5, GHQ_w6, GHQ_w7, GHQ_w8, GHQ_w9, GHQ_w10,
2336 :     GHQ_w11};
2337 :    
2338 :     /**
2339 :     * Evaluate the conditional deviance for a given set of fixed effects.
2340 :     *
2341 :     * @param GS Pointer to a GlmerStruct
2342 :     * @param fixed value of the fixed effects
2343 :     *
2344 :     * @return conditional deviance
2345 :     */
2346 :     static double
2347 :     fixed_effects_deviance(GlmerStruct GS, const double fixed[])
2348 :     {
2349 :     SEXP devs;
2350 :     int i, ione = 1;
2351 :     double ans, one = 1, zero = 0;
2352 :    
2353 :     F77_CALL(dgemv)("N", &(GS->n), &(GS->p), &one,
2354 :     GS->X, &(GS->n), fixed,
2355 :     &ione, &zero, REAL(GS->eta), &ione);
2356 :     /* add in random effects and offset */
2357 :     vecIncrement(REAL(GS->eta), GS->off, GS->n);
2358 :     eval_check_store(GS->linkinv, GS->rho, GS->mu);
2359 :     devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));
2360 :     for (i = 0, ans = 0; i < GS->n; i++) ans += REAL(devs)[i];
2361 :     UNPROTECT(1);
2362 :     return ans;
2363 :     }
2364 :    
2365 :    
2366 :     /**
2367 :     * Evaluate the quadratic form (x-mn)'A'A(x-mn) from the multivariate
2368 :     * normal kernel.
2369 :     *
2370 :     * @param n dimension of random variate
2371 :     * @param mn mean vector
2372 :     * @param a upper Cholesky factor of the precision matrix
2373 :     * @param lda leading dimension of A
2374 :     * @param x vector at which to evaluate the kernel
2375 :     *
2376 :     * @return value of the normal kernel
2377 :     */
2378 :     static double
2379 :     normal_kernel(int n, const double mn[],
2380 :     const double a[], int lda, const double x[])
2381 :     {
2382 :     int i, ione = 1;
2383 :     double *tmp = Calloc(n, double), ans;
2384 :    
2385 :     for (i = 0; i < n; i++) tmp[i] = x[i] - mn[i];
2386 :     F77_CALL(dtrmv)("U", "N", "N", &n, a, &lda, tmp, &ione);
2387 :     for (i = 0, ans = 0; i < n; i++) ans += tmp[i] * tmp[i];
2388 :     Free(tmp);
2389 :     return ans;
2390 :     }
2391 :    
2392 :     /* FIXME: Separate the calculation of the offset random effects from
2393 :     * the calculation of the deviance contributions. */
2394 :    
2395 :     /**
2396 :     * Determine the deviance components associated with each of the
2397 :     * levels of a grouping factor at the conditional modes or a value
2398 :     * offset from the conditional modes by delb.
2399 :     *
2400 :     * @param GS pointer to a GlmerStruct
2401 :     * @param b conditional modes of the random effects
2402 :     * @param Gp group pointers
2403 :     * @param nc number of columns in the model matrix for the kth
2404 :     * grouping factor
2405 :     * @param k index (0-based) of the grouping factor
2406 :     * @param delb vector of length nc giving the changes in the
2407 :     * orthonormalized random effects
2408 :     * @param OmgFac Cholesky factor of the inverse of the penalty matrix
2409 :     * for this grouping factor
2410 :     * @param bVfac 3-dimensional array holding the factors of the
2411 :     * conditional variance-covariance matrix of the random effects
2412 :     FIXME: This is wrong. It is bVar[[i]] not bVfac that is being passed.
2413 :     This only affects the AGQ method.
2414 :     * @param devcmp array to hold the deviance components
2415 :     *
2416 :     * @return devcmp
2417 :     */
2418 :     static double*
2419 :     rel_dev_1(GlmerStruct GS, const double b[], int nlev, int nc, int k,
2420 :     const double delb[], const double OmgFac[],
2421 :     const double bVfac[], double devcmp[])
2422 :     {
2423 :     SEXP devs;
2424 :     int *fv = INTEGER(VECTOR_ELT(GET_SLOT(GS->mer, Matrix_flistSym), k)),
2425 :     i, j;
2426 :     double *bcp = (double *) NULL;
2427 :    
2428 :     AZERO(devcmp, nlev);
2429 :     if (delb) {
2430 :     int ione = 1, ntot = nlev * nc;
2431 :     double sumsq = 0;
2432 :     /* copy the contents of b */
2433 :     bcp = Memcpy(Calloc(ntot, double), b, ntot);
2434 :     if (nc == 1) {
2435 :     sumsq = delb[0] * delb[0];
2436 :     for (i = 0; i < nlev; i++) b[i] += delb[0] * bVfac[i];
2437 :     } else {
2438 :     int ncsq = nc * nc;
2439 :     double *tmp = Calloc(nc, double);
2440 :     for (i = 0; i < nlev; i++) {
2441 :     Memcpy(tmp, delb, nc);
2442 :     F77_CALL(dtrmv)("U", "N", "N", &nc, &(bVfac[i * ncsq]),
2443 :     &nc, tmp, &ione);
2444 :     for (j = 0; j < nc; j++) b[i + j * nc] = tmp[j];
2445 :     }
2446 :     /* sum of squares of delb */
2447 :     for (j = 0; j < nc; j++) sumsq += delb[j] * delb[j];
2448 :     }
2449 :     for (i = 0; i < nlev; i++) devcmp[i] = -sumsq;
2450 :     }
2451 :     internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
2452 :     eval_check_store(GS->linkinv, GS->rho, GS->mu);
2453 :     devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));
2454 :     for (i = 0; i < GS->n; i++)
2455 :     devcmp[fv[i] - 1] += REAL(devs)[i];
2456 :     UNPROTECT(1);
2457 :     if (nc == 1) {
2458 :     for (i = 0; i < nlev; i++) {
2459 :     double tmp = *OmgFac * b[i];
2460 :     devcmp[i] += tmp * tmp;
2461 :     }
2462 :     } else {
2463 :     double *tmp = Calloc(nc, double);
2464 :     int ione = 1;
2465 :    
2466 :     for (i = 0; i < nlev; i++) {
2467 :     for (j = 0; j < nc; j++) tmp[j] = b[i + j * nlev];
2468 :     F77_CALL(dtrmv)("U", "N", "N", &nc, OmgFac, &nc,
2469 :     tmp, &ione);
2470 :     for (j = 0; j < nc; j++)
2471 :     devcmp[i] += tmp[j] * tmp[j];
2472 :     }
2473 :     }
2474 :     if (delb) {
2475 :     Memcpy(b, bcp, ntot);
2476 :     Free(bcp);
2477 :     }
2478 :     return devcmp;
2479 :     }
2480 :    
2481 :    
2482 :     /**
2483 :     * Create a Markov Chain Monte Carlo sample from a fitted generalized
2484 :     * linear mixed model
2485 :     *
2486 :     * @param GSpt External pointer to a GlmerStruct
2487 :     * @param b Conditional modes of the random effects at the parameter
2488 :     * estimates
2489 :     * @param fixed Estimates of the fixed effects
2490 :     * @param varc Estimates of the variance components
2491 :     * @param savebp Logical indicator of whether or not to save the
2492 :     * random effects in the MCMC sample
2493 :     * @param nsampp Integer value of the number of samples to generate
2494 :     *
2495 :     * @return
2496 :     */
2497 :     SEXP
2498 :     glmer_MCMCsamp(SEXP GSpt, SEXP b, SEXP fixed, SEXP varc,
2499 :     SEXP savebp, SEXP nsampp)
2500 :     {
2501 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSpt);
2502 :     int i, j, nf = LENGTH(b), nsamp = asInteger(nsampp),
2503 :     p = LENGTH(fixed), q = LENGTH(varc),
2504 :     saveb = asLogical(savebp);
2505 :     int *mcpar, nc = p + q;
2506 :     SEXP ans, mcparSym = install("mcpar");
2507 :    
2508 :     if (nsamp <= 0) nsamp = 1;
2509 :     nc = p + q;
2510 :     if (saveb)
2511 :     for (i = 0; i < nf; i++) {
2512 :     int *dims = INTEGER(getAttrib(VECTOR_ELT(b, i),
2513 :     R_DimSymbol));
2514 :     nc += dims[0] * dims[1];
2515 :     }
2516 :     ans = PROTECT(allocMatrix(REALSXP, nsamp, nc));
2517 :     GetRNGstate();
2518 :     for (i = 0; i < nsamp; i++) {
2519 :     internal_glmer_fixef_update(GS, b, REAL(fixed));
2520 :     internal_bhat(GS, REAL(fixed), REAL(varc));
2521 :     internal_glmer_ranef_update(GS, b);
2522 :     internal_Omega_update(GS->mer, b);
2523 :     internal_mer_coef(GS->mer, 2, REAL(varc));
2524 :     for (j = 0; j < p; j++)
2525 :     REAL(ans)[i + j * nsamp] = REAL(fixed)[j];
2526 :     for (j = 0; j < q; j++)
2527 :     REAL(ans)[i + (j + p) * nsamp] = REAL(varc)[j];
2528 :     if (saveb) {
2529 :     int base = p + q, k;
2530 :     for (k = 0; k < nf; k++) {
2531 :     SEXP bk = VECTOR_ELT(b, k);
2532 :     int *dims = INTEGER(getAttrib(bk, R_DimSymbol));
2533 :     int klen = dims[0] * dims[1];
2534 :    
2535 :     for (j = 0; j < klen; j++)
2536 :     REAL(ans)[i + (j + base) * nsamp] = REAL(bk)[j];
2537 :     base += klen;
2538 :     }
2539 :     }
2540 :     }
2541 :     PutRNGstate();
2542 :     UNPROTECT(1);
2543 :     /* set (S3) class and mcpar attribute */
2544 :     setAttrib(ans, R_ClassSymbol, mkString("mcmc"));
2545 :     setAttrib(ans, mcparSym, allocVector(INTSXP, 3));
2546 :     mcpar = INTEGER(getAttrib(ans, mcparSym));
2547 :     mcpar[0] = mcpar[2] = 1;
2548 :     mcpar[1] = nsamp;
2549 :    
2550 :     return ans;
2551 :     }
2552 :    
2553 :     /**
2554 :     * Determine the conditional modes and the conditional variance of the
2555 :     * fixed effects given the data and the current random effects.
2556 :     * Create a Metropolis-Hasting proposal step from the multivariate
2557 :     * normal density, determine the acceptance probability and return the
2558 :     * current value or the proposed value.
2559 :     *
2560 :     * @param GSp pointer to a GlmerStruct
2561 :     * @param b list of random effects
2562 :     * @param fixed current value of the fixed effects
2563 :     *
2564 :     * @return updated value of the fixed effects
2565 :     */
2566 :     SEXP glmer_fixed_update(SEXP GSp, SEXP b, SEXP fixed)
2567 :     {
2568 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
2569 :    
2570 :     if (!isReal(fixed) || LENGTH(fixed) != GS->p)
2571 :     error(_("%s must be a %s of length %d"), "fixed",
2572 :     "numeric vector", GS->p);
2573 :     GetRNGstate();
2574 :     internal_glmer_fixef_update(GS, b, REAL(fixed));
2575 :     PutRNGstate();
2576 :     return fixed;
2577 :     }
2578 :    
2579 :     SEXP glmer_bhat(SEXP pars, SEXP GSp)
2580 :     {
2581 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
2582 :    
2583 :     if (!isReal(pars) || LENGTH(pars) != GS->npar)
2584 :     error(_("`%s' must be a numeric vector of length %d"),
2585 :     "pars", GS->npar);
2586 :     if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))
2587 :     warning(_("internal_bhat did not converge"));
2588 :     return R_NilValue;
2589 :     }
2590 :    
2591 :    
2592 :    
2593 :     /**
2594 :     * Determine the conditional modes and the conditional variance of the
2595 :     * fixed effects given the data and the current random effects.
2596 :     * Create a Metropolis-Hasting proposal step from the multivariate
2597 :     * normal density, determine the acceptance probability and, if the
2598 :     * step is to be accepted, overwrite the contents of fixed with the
2599 :     * new contents.
2600 :     *
2601 :     * @param GS a GlmerStruct
2602 :     * @param b list of random effects
2603 :     * @param fixed current value of the fixed effects
2604 :     *
2605 :     * @return updated value of the fixed effects
2606 :     */
2607 :     static double *
2608 :     internal_glmer_fixef_update(GlmerStruct GS, SEXP b,
2609 :     double fixed[])
2610 :     {
2611 :     SEXP dmu_deta, var;
2612 :     int i, ione = 1, it, j, lwork = -1;
2613 :     double *ans = Calloc(GS->p, double), /* proposal point */
2614 :     *md = Calloc(GS->p, double), /* conditional modes */
2615 :     *w = Calloc(GS->n, double), *work,
2616 :     *wtd = Calloc(GS->n * GS->p, double),
2617 :     *z = Calloc(GS->n, double),
2618 :     crit, devr, one = 1, tmp, zero = 0;
2619 :    
2620 :     if (!isNewList(b) || LENGTH(b) != GS->nf)
2621 :     error(_("%s must be a %s of length %d"), "b", "list", GS->nf);
2622 :     for (i = 0; i < GS->nf; i++) {
2623 :     SEXP bi = VECTOR_ELT(b, i);
2624 :     if (!isReal(bi) || !isMatrix(bi))
2625 :     error(_("b[[%d]] must be a numeric matrix"), i);
2626 :     }
2627 :     AZERO(z, GS->n); /* -Wall */
2628 :     Memcpy(md, fixed, GS->p);
2629 :     /* calculate optimal size of work array */
2630 :     F77_CALL(dgels)("N", &(GS->n), &(GS->p), &ione, wtd, &(GS->n),
2631 :     z, &(GS->n), &tmp, &lwork, &j);
2632 :     if (j) /* shouldn't happen */
2633 :     error(_("%s returned error code %d"), "dgels", j);
2634 :     lwork = (int) tmp;
2635 :     work = Calloc(lwork, double);
2636 :    
2637 :     AZERO(GS->off, GS->n); /* fitted values from random effects */
2638 :     /* fitted_ranef(GET_SLOT(GS->mer, Matrix_flistSym), GS->unwtd, b, */
2639 :     /* INTEGER(GET_SLOT(GS->mer, Matrix_ncSym)), GS->off); */
2640 :     for (i = 0; i < GS->n; i++)
2641 :     (GS->etaold)[i] = ((GS->off)[i] += (GS->offset)[i]);
2642 :    
2643 :     for (it = 0, crit = GS->tol + 1;
2644 :     it < GS->maxiter && crit > GS->tol; it++) {
2645 :     /* fitted values from current beta */
2646 :     F77_CALL(dgemv)("N", &(GS->n), &(GS->p), &one,
2647 :     GS->X, &(GS->n), md,
2648 :     &ione, &zero, REAL(GS->eta), &ione);
2649 :     /* add in random effects and offset */
2650 :     vecIncrement(REAL(GS->eta), (GS->off), GS->n);
2651 :     /* check for convergence */
2652 :     crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);
2653 :     /* obtain mu, dmu_deta, var */
2654 :     eval_check_store(GS->linkinv, GS->rho, GS->mu);
2655 :     dmu_deta = PROTECT(eval_check(GS->mu_eta, GS->rho,
2656 :     REALSXP, GS->n));
2657 :     var = PROTECT(eval_check(GS->var, GS->rho, REALSXP, GS->n));
2658 :     /* calculate weights and working residual */
2659 :     for (i = 0; i < GS->n; i++) {
2660 :     w[i] = GS->wts[i] * REAL(dmu_deta)[i]/sqrt(REAL(var)[i]);
2661 :     z[i] = w[i] * (REAL(GS->eta)[i] - (GS->off)[i] +
2662 :     ((GS->y)[i] - REAL(GS->mu)[i]) /
2663 :     REAL(dmu_deta)[i]);
2664 :     }
2665 :     UNPROTECT(2);
2666 :     /* weighted copy of the model matrix */
2667 :     for (j = 0; j < GS->p; j++)
2668 :     for (i = 0; i < GS->n; i++)
2669 :     wtd[i + j * GS->n] = GS->X[i + j * GS->n] * w[i];
2670 :     /* weighted least squares solution */
2671 :     F77_CALL(dgels)("N", &(GS->n), &(GS->p), &ione, wtd, &(GS->n),
2672 :     z, &(GS->n), work, &lwork, &j);
2673 :     if (j) error(_("%s returned error code %d"), "dgels", j);
2674 :     Memcpy(md, z, GS->p);
2675 :     }
2676 :     /* wtd contains the Cholesky factor of
2677 :     * the precision matrix */
2678 :     devr = normal_kernel(GS->p, md, wtd, GS->n, fixed);
2679 :     devr -= fixed_effects_deviance(GS, fixed);
2680 :     for (i = 0; i < GS->p; i++) {
2681 :     double var = norm_rand();
2682 :     ans[i] = var;
2683 :     devr -= var * var;
2684 :     }
2685 :     F77_CALL(dtrsv)("U", "N", "N", &(GS->p), wtd, &(GS->n), ans, &ione);
2686 :     for (i = 0; i < GS->p; i++) ans[i] += md[i];
2687 :     devr += fixed_effects_deviance(GS, ans);
2688 :     crit = exp(-0.5 * devr); /* acceptance probability */
2689 :     tmp = unif_rand();
2690 :     if (asLogical(Matrix_getElement(GS->cv, "msVerbose"))) {
2691 :     Rprintf("%5.3f: ", crit);
2692 :     for (j = 0; j < GS->p; j++) Rprintf("%#10g ", ans[j]);
2693 :     Rprintf("\n");
2694 :     }
2695 :     if (tmp < crit) Memcpy(fixed, ans, GS->p);
2696 :     Free(ans); Free(md); Free(w);
2697 :     Free(work); Free(wtd); Free(z);
2698 :     return fixed;
2699 :     }
2700 :    
2701 :     static void
2702 :     internal_glmer_ranef_update(GlmerStruct GS, SEXP b)
2703 :     {
2704 :     SEXP bhat, bprop = PROTECT(duplicate(b)),
2705 :     bVar = GET_SLOT(GS->mer, Matrix_bVarSym),
2706 :     Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);
2707 :     int i, ione = 1, j, k;
2708 :     double *b = REAL(GET_SLOT(GS->mer, Matrix_ranefSym);
2709 :     double devr, one = 1;
2710 :    
2711 :     bhat = R_NilValue;
2712 :     /* subtract deviance at b */
2713 :     devr = -random_effects_deviance(GS, b);
2714 :     for (i = 0; i < GS->nf; i++) {
2715 :     SEXP Bi = VECTOR_ELT(b, i);
2716 :     int *dims = INTEGER(getAttrib(Bi, R_DimSymbol));
2717 :     int nlev = dims[0], nci = dims[1];
2718 :     int ncisqr = nci * nci, ntot = nlev * nci;
2719 :     double *bcopy = Calloc(ntot, double),
2720 :     *bi = REAL(Bi),
2721 :     *bhati = REAL(VECTOR_ELT(bhat, i)),
2722 :     *bpropi = REAL(VECTOR_ELT(bprop, i)),
2723 :     *bVari = REAL(VECTOR_ELT(bVar, i)),
2724 :     *chol = Calloc(ncisqr, double),
2725 :     *delta = Calloc(nci, double),
2726 :     *omgfac = Memcpy(Calloc(ncisqr, double),
2727 :     REAL(VECTOR_ELT(Omega, i)),
2728 :     ncisqr);
2729 :    
2730 :     /* subtract quadratic form in */
2731 :     /* Omega[[i]] at b */
2732 :     F77_CALL(dpotrf)("U", &nci, omgfac, &nci, &j);
2733 :     if (j)
2734 :     error(_("Leading %d minor of Omega[[%d]] not positive definite"),
2735 :     j, i + 1);
2736 :     Memcpy(bcopy, bi, ntot);
2737 :     F77_CALL(dtrmm)("R", "U", "T", "N", &nlev, &nci, &one,
2738 :     omgfac, &nci, bcopy, &nlev);
2739 :     for (k = 0; k < ntot; k++) devr -= bcopy[k] * bcopy[k];
2740 :     /* form bprop and proposal density */
2741 :     for (k = 0; k < nlev; k++) {
2742 :     /* proposal density at b */
2743 :     for (j = 0; j < nci; j++)
2744 :     delta[j] = bi[k + j * nlev] - bhati[k + j * nlev];
2745 :     Memcpy(chol, &(bVari[k * ncisqr]), ncisqr);
2746 :     F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);
2747 :     if (j)
2748 :     error(_("Leading %d minor of bVar[[%d]][,,%d] not positive definite"),
2749 :     j, i + 1, k + 1);
2750 :     F77_CALL(dtrsv)("U", "T", "N", &nci, chol, &nci,
2751 :     delta, &ione);
2752 :     for (j = 0; j < nci; j++) {
2753 :     double nrm = norm_rand(); /* proposal deviate */
2754 :     devr += delta[j] * delta[j] - nrm * nrm;
2755 :     delta[j] = nrm;
2756 :     }
2757 :     /* scale by Cholesky inverse */
2758 :     F77_CALL(dtrmv)("U", "T", "N", &nci, chol, &nci,
2759 :     delta, &ione);
2760 :     /* add mean */
2761 :     for (j = 0; j < nci; j++)
2762 :     bpropi[k + j * nlev] = bhati[k + j * nlev] + delta[j];
2763 :     }
2764 :     /* add quadratic form in */
2765 :     /* Omega[[i]] at bprop */
2766 :     Memcpy(bcopy, bpropi, ntot);
2767 :     F77_CALL(dtrmm)("R", "U", "T", "N", &nlev, &nci, &one,
2768 :     omgfac, &nci, bcopy, &nlev);
2769 :     for (k = 0; k < ntot; k++) devr += bcopy[k] * bcopy[k];
2770 :    
2771 :     Free(bcopy); Free(chol); Free(delta); Free(omgfac);
2772 :     }
2773 :     /* add deviance at bprop */
2774 :     /* devr += random_effects_deviance(GS, bprop); */
2775 :    
2776 :     if (unif_rand() < exp(-0.5 * devr))
2777 :     for (i = 0; i < GS->nf; i++) { /* copy each face of b */
2778 :     SEXP Bi = VECTOR_ELT(b, i);
2779 :     int *dims = INTEGER(getAttrib(Bi, R_DimSymbol));
2780 :    
2781 :     Memcpy(REAL(Bi), REAL(VECTOR_ELT(bprop, i)),
2782 :     dims[0] * dims[1]);
2783 :     }
2784 :    
2785 :     if (asLogical(Matrix_getElement(GS->cv, "msVerbose"))) {
2786 :     double *b0 = REAL(VECTOR_ELT(bprop, 0));
2787 :     Rprintf("%5.3f:", exp(-0.5 * devr));
2788 :     for (k = 0; k < 5; k++) Rprintf("%#10g ", b0[k]);
2789 :     Rprintf("\n");
2790 :     }
2791 :    
2792 :     UNPROTECT(2);
2793 :     }
2794 :    
2795 :    
2796 :     /**
2797 :     * Compute the approximation to the deviance using adaptive
2798 :     * Gauss-Hermite quadrature (AGQ). When nAGQ == 1 this is the Laplace
2799 :     * approximation.
2800 :     *
2801 :     * @param pars pointer to a numeric vector of parameters
2802 :     * @param GSp pointer to a GlmerStruct object
2803 :     * @param nAGQp pointer to a scalar integer representing the number of
2804 :     * points in AGQ to use
2805 :     *
2806 :     * @return the approximation to the deviance as computed using AGQ
2807 :     */
2808 :     SEXP glmer_devAGQ(SEXP pars, SEXP GSp, SEXP nAGQp)
2809 :     {
2810 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
2811 :     SEXP Omega = GET_SLOT(GS->mer, Matrix_OmegaSym),
2812 :     bVar = GET_SLOT(GS->mer, Matrix_bVarSym);
2813 :     int i, j, k, nAGQ = asInteger(nAGQp);
2814 :     int n2 = (nAGQ + 1)/2,
2815 :     *Gp = INTEGER(GET_SLOT(GS->mer, Matrix_GpSym)),
2816 :     *nc = INTEGER(GET_SLOT(GS->mer, Matrix_ncSym));
2817 :     double *f0, LaplaceDev = 0, AGQadjst = 0,
2818 :     *bhat = REAL(GET_SLOT(GS->mer, Matrix_ranefSym));
2819 :    
2820 :     if (!isReal(pars) || LENGTH(pars) != GS->npar)
2821 :     error(_("`%s' must be a numeric vector of length %d"),
2822 :     "pars", GS->npar);
2823 :     if (GS->nf > 1 && nAGQ > 1) {
2824 :     warning(_("AGQ not available for multiple grouping factors - using Laplace"));
2825 :     nAGQ = 1;
2826 :     }
2827 :     if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))
2828 :     return ScalarReal(DBL_MAX);
2829 :    
2830 :     for (i = 0; i < GS->nf; i++) {
2831 :     int nci = nc[i];
2832 :     int ncip1 = nci + 1, ncisqr = nci * nci,
2833 :     nlev = (Gp[i + 1] - Gp[i])/nci;
2834 :     double *omgf = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)), Matrix_xSym)),
2835 :     *bVi = Memcpy(Calloc(ncisqr * nlev, double),
2836 :     REAL(VECTOR_ELT(bVar, i)), ncisqr * nlev);
2837 :    
2838 :     for (j = 0; j < nci; j++) { /* nlev * logDet(Omega_i) */
2839 :     LaplaceDev += 2 * nlev * log(omgf[j * ncip1]);
2840 :     }
2841 :     for (k = 0; k < nlev; k++) {
2842 :     double *bVik = bVi + k * ncisqr;
2843 :     F77_CALL(dpotrf)("U", &nci, bVik, &nci, &j);
2844 :     if (j)
2845 :     error(_("Leading %d minor of bVar[[%d]][,,%d] not positive definite"),
2846 :     j, i + 1, k + 1);
2847 :     for (j = 0; j < nci; j++) LaplaceDev -= 2 * log(bVik[j * ncip1]);
2848 :     }
2849 :    
2850 :     f0 = Calloc(nlev, double);
2851 :     rel_dev_1(GS, bhat, nlev, nci, i, (double *) NULL,
2852 :     omgf, bVi, f0);
2853 :     for (k = 0; k < nlev; k++) LaplaceDev += f0[k];
2854 :     if (nAGQ > 1) {
2855 :     double *fx = Calloc(nlev, double),
2856 :     *rellik = Calloc(nlev, double),
2857 :     *delb = Calloc(nci, double);
2858 :    
2859 :     if (nci > 1) error(_("code not yet written"));
2860 :     AZERO(rellik, nlev); /* zero accumulator */
2861 :     for (k = 0; k < n2; k++) {
2862 :     delb[0] = GHQ_x[nAGQ][k];
2863 :     if (delb[0]) {
2864 :     rel_dev_1(GS, bhat, nlev, nci, i, delb,
2865 :     omgf, bVi, fx);
2866 :     for (j = 0; j < nlev; j++) {
2867 :     rellik[j] += GHQ_w[nAGQ][k] *
2868 :     exp(-(fx[j] - f0[j])/2);
2869 :     }
2870 :     delb[0] *= -1;
2871 :     rel_dev_1(GS, bhat, nlev, nci, i, delb,
2872 :     omgf, bVi, fx);
2873 :     for (j = 0; j < nlev; j++) {
2874 :     rellik[j] += GHQ_w[nAGQ][k] *
2875 :     exp(-(fx[j] - f0[j])/2);
2876 :     }
2877 :     } else {
2878 :     for (j = 0; j < nlev; j++)
2879 :     rellik[j] += GHQ_w[nAGQ][k];
2880 :     }
2881 :     }
2882 :     for (j = 0; j < nlev; j++)
2883 :     AGQadjst -= 2 * log(rellik[j]);
2884 :     Free(fx); Free(rellik);
2885 :     }
2886 :     Free(f0); Free(bVi);
2887 :     }
2888 :     return ScalarReal(LaplaceDev + AGQadjst);
2889 :     }
2890 :    
2891 :    
2892 :     /**
2893 :     * Determine the conditional modes and the conditional variance of the
2894 :     * random effects given the data and the current fixed effects and
2895 :     * variance components.
2896 :     *
2897 :     * Create a Metropolis-Hasting proposal step from a multivariate
2898 :     * normal density centered at bhat with variance-covariance matrix
2899 :     * from bVar, determine the acceptance probability and return the
2900 :     * current value or the proposed value.
2901 :     *
2902 :     * @param GSp pointer to a GlmerStruct
2903 :     * @param fixed pointer to a numeric vector of the fixed effects
2904 :     * @param varc pointer to a numeric vector of the variance components
2905 :     * @param varc pointer to current values of b
2906 :     *
2907 :     * @return updated b from the Metropolis-Hastings step
2908 :     */
2909 :     SEXP glmer_ranef_update(SEXP GSp, SEXP fixed, SEXP varc, SEXP b)
2910 :     {
2911 :     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
2912 :     int nvarc = GS->npar - GS->p;
2913 :    
2914 :     if (!isReal(fixed) || LENGTH(fixed) != GS->p)
2915 :     error(_("`%s' must be a numeric vector of length %d"),
2916 :     "fixed", GS->p);
2917 :     if (INTEGER(GET_SLOT(GS->mer, Matrix_ncSym))[GS->nf] > 0)
2918 :     error(_("the mer object must be set to skip fixed effects"));
2919 :     if (!isReal(varc) || LENGTH(varc) != nvarc)
2920 :     error(_("`%s' must be a numeric vector of length %d"),
2921 :     "varc", nvarc);
2922 :    
2923 :     GetRNGstate();
2924 :     /* Don't check for convergence failure after internal_bhat.
2925 :     * It is determining the mean of the proposal density and
2926 :     * does not need exact convergence. */
2927 :     internal_bhat(GS, REAL(fixed), REAL(varc));
2928 :     internal_glmer_ranef_update(GS, b);
2929 :     PutRNGstate();
2930 :    
2931 :     UNPROTECT(1);
2932 :     return b;
2933 :     }
2934 :    
2935 :    
2936 :     #endif

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge