SCM

SCM Repository

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

Diff of /pkg/src/lmer.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1161, Fri Jan 13 23:46:29 2006 UTC revision 1162, Fri Jan 13 23:47:38 2006 UTC
# Line 47  Line 47 
47      return ans;      return ans;
48  }  }
49    
 /* These should be moved to other source files such as chm_common.c */  
 static cholmod_dense  
 *numeric_as_chm_dense(double *v, int n)  
 {  
     cholmod_dense *ans = Calloc(1, cholmod_dense);  
   
     ans->d = ans->nzmax = ans->nrow = n;  
     ans->ncol = 1;  
     ans->x = (void *) v;  
     ans->xtype = CHOLMOD_REAL;  
     ans->dtype = CHOLMOD_DOUBLE;  
     return ans;  
 }  
   
 static  
 SEXP alloc_dgeMatrix(int m, int n, SEXP rownms, SEXP colnms)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dgeMatrix"))), dn;  
     int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  
   
     dims[0] = m; dims[1] = n;  
     ALLOC_SLOT(ans, Matrix_xSym, REALSXP, m * n);  
     dn = GET_SLOT(ans, Matrix_DimNamesSym);  
     SET_VECTOR_ELT(dn, 0, duplicate(rownms));  
     SET_VECTOR_ELT(dn, 1, duplicate(colnms));  
     UNPROTECT(1);  
     return ans;  
 }  
   
 static  
 SEXP alloc_dpoMatrix(int n, char *uplo, SEXP rownms, SEXP colnms)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dpoMatrix"))), dn;  
     int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  
   
     dims[0] = dims[1] = n;  
     ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * n);  
     SET_SLOT(ans, Matrix_uploSym, mkString(uplo));  
     dn = GET_SLOT(ans, Matrix_DimNamesSym);  
     SET_VECTOR_ELT(dn, 0, duplicate(rownms));  
     SET_VECTOR_ELT(dn, 1, duplicate(colnms));  
     UNPROTECT(1);  
     return ans;  
 }  
   
 static  
 SEXP alloc_dtrMatrix(int n, char *uplo, char *diag, SEXP rownms, SEXP colnms)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dtrMatrix"))), dn;  
     int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  
   
     dims[0] = dims[1] = n;  
     ALLOC_SLOT(ans, Matrix_xSym, REALSXP, n * n);  
     SET_SLOT(ans, Matrix_uploSym, mkString(uplo));  
     SET_SLOT(ans, Matrix_diagSym, mkString(diag));  
     dn = GET_SLOT(ans, Matrix_DimNamesSym);  
     SET_VECTOR_ELT(dn, 0, duplicate(rownms));  
     SET_VECTOR_ELT(dn, 1, duplicate(colnms));  
     UNPROTECT(1);  
     return ans;  
 }  
   
 static  
 SEXP alloc_dsCMatrix(int n, int nz, char *uplo, SEXP rownms, SEXP colnms)  
 {  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("dsCMatrix"))), dn;  
     int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2));  
   
     dims[0] = dims[1] = n;  
     ALLOC_SLOT(ans, Matrix_xSym, REALSXP, nz);  
     ALLOC_SLOT(ans, Matrix_iSym, INTSXP, nz);  
     ALLOC_SLOT(ans, Matrix_pSym, INTSXP, n + 1);  
     SET_SLOT(ans, Matrix_uploSym, mkString(uplo));  
     dn = GET_SLOT(ans, Matrix_DimNamesSym);  
     SET_VECTOR_ELT(dn, 0, duplicate(rownms));  
     SET_VECTOR_ELT(dn, 1, duplicate(colnms));  
     UNPROTECT(1);  
     return ans;  
 }  
50    
51  /* Internally used utilities */  /* Internally used utilities */
52    
# Line 193  Line 114 
114      }      }
115  }  }
116    
117    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  /**  /**
139   * Evaluate the quadratic form in b defined by Omega   * Evaluate the quadratic form in b defined by Omega
140   *   *
# Line 330  Line 272 
272      SEXP rho;        /* environment in which to evaluate the calls */      SEXP rho;        /* environment in which to evaluate the calls */
273      SEXP eta;        /* linear predictor */      SEXP eta;        /* linear predictor */
274      SEXP mu;         /* mean vector */      SEXP mu;         /* mean vector */
     SEXP unwtd;      /* unweighted model matrices */  
     SEXP wtd;        /* weighted model matrices */  
275      SEXP linkinv;    /* expression for inverse link evaluation */      SEXP linkinv;    /* expression for inverse link evaluation */
276      SEXP mu_eta;     /* expression for dmu/deta evaluation */      SEXP mu_eta;     /* expression for dmu/deta evaluation */
277      SEXP LMEopt;     /* expression for LME optimization */      SEXP LMEopt;     /* expression for LME optimization */
278      SEXP dev_resids; /* expression for deviance residuals */      SEXP dev_resids; /* expression for deviance residuals */
279      SEXP var;        /* expression for variance evaluation */      SEXP var;        /* expression for variance evaluation */
     double *off;     /* offset for random effects only */  
280      double *offset;  /* offset for GLM */      double *offset;  /* offset for GLM */
281      double *wts;     /* prior weights for GLM */      double *wts;     /* prior weights for GLM */
     double *w;       /* vector of weights */  
     double *X;       /* copy of fixed-effects model matrix */  
282      double *y;       /* copy of response vector */      double *y;       /* copy of response vector */
     double *z;       /* working residual */  
     double *Zt;      /* copy of x slot in transpose of random-effects model matrix */  
283      double *etaold;  /* previous value of eta for evaluating  */      double *etaold;  /* previous value of eta for evaluating  */
284      int n;           /* length of the response vector */      int n;           /* length of the response vector */
285      int p;           /* length of fixed effects vector */      int p;           /* length of fixed effects vector */
# Line 356  Line 291 
291      double tol;      /* convergence tolerance for IRLS iterations */      double tol;      /* convergence tolerance for IRLS iterations */
292  } glmer_struct, *GlmerStruct;  } glmer_struct, *GlmerStruct;
293    
 static void  
 internal_mer_bVar(SEXP x)  
 {  
     int q = LENGTH(GET_SLOT(x, Matrix_rZySym));  
     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));  
     cholmod_sparse *eye = cholmod_speye(q, q, CHOLMOD_REAL, &c), *Linv;  
     int *Perm = (int *)(L->Perm), *iperm = Calloc(q, int), i;  
                                 /* create the inverse permutation */  
     for (i = 0; i < q; i++) iperm[Perm[i]] = i;  
                                 /* apply iperm to the identity matrix */  
     for (i = 0; i < q; i++) ((int*)(eye->i))[i] = iperm[i];  
                                 /* Create Linv */  
     Linv = cholmod_spsolve(CHOLMOD_L, L, eye, &c);  
     cholmod_free_sparse(&eye, &c);  
     Linv_to_bVar(Linv, INTEGER(GET_SLOT(x, Matrix_GpSym)),  
                  INTEGER(GET_SLOT(x, Matrix_ncSym)),  
                  GET_SLOT(x, Matrix_bVarSym), "U");  
     cholmod_free_sparse(&Linv, &c);  
     Free(L);  
 }  
   
294  /**  /**
295   * Increment the fitted values from an mer object using the (possibly   * Calculate fitted values for the current fixed and random effects.
  * externally stored) contents of the model matrices.  
296   *   *
297   * @param x pointer to an mer object   * @param x pointer to an mer object
  * @param X contents of the X model matrix or (double *) NULL  
  * @param Zt values of non-zeros in the Zt model matrix or (double *) NULL  
298   * @param initial initial values (i.e. an offset) or (double *) NULL   * @param initial initial values (i.e. an offset) or (double *) NULL
299   * @param val array to hold the fitted values   * @param val array to hold the fitted values
300   *   *
301   * @return pointer to a numeric array of fitted values   * @return pointer to a numeric array of fitted values
302   */   */
303  static double *  static double *
304  internal_mer_fitted(SEXP x, const double X[], const double Ztx[],  internal_mer_fitted(SEXP x, const double initial[], double val[])
                     const double initial[], double val[])  
305  {  {
306      int ione = 1, n = LENGTH(GET_SLOT(x, Matrix_ySym));      SEXP fixef = GET_SLOT(x, Matrix_fixefSym),
307      double one[] = {1,0};          ranef = GET_SLOT(x, Matrix_ranefSym);
308        int ione = 1, n = LENGTH(GET_SLOT(x, Matrix_ySym)), p = LENGTH(fixef);
309      mer_secondary(x);      double *X = REAL(GET_SLOT(x, Matrix_XSym)), one[] = {1,0};
     if (initial) Memcpy(val, initial, n);  
     else AZERO(val, n);  
     if (X) {  
         SEXP fixef = GET_SLOT(x, Matrix_fixefSym);  
         int p = LENGTH(fixef);  
   
         F77_CALL(dgemv)("N", &n, &p, one, X, &n, REAL(fixef),  
                         &ione, one, val, &ione);  
     }  
     if (Ztx) {  
         SEXP ranef = GET_SLOT(x, Matrix_ranefSym);  
310          cholmod_sparse *Zt = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtSym));          cholmod_sparse *Zt = as_cholmod_sparse(GET_SLOT(x, Matrix_ZtSym));
311          cholmod_dense *chv = numeric_as_chm_dense(val, n),          cholmod_dense *chv = numeric_as_chm_dense(val, n),
312              *chb = as_cholmod_dense(ranef);              *chb = as_cholmod_dense(ranef);
         cholmod_sparse *Ztcp = cholmod_copy_sparse(Zt, &c);  
313    
314          if (Ztx != (double*)(Zt->x)) Memcpy((double*)(Ztcp->x), Ztx, cholmod_nnz(Zt, &c));      mer_secondary(x);
315          Free(Zt);      if (initial) Memcpy(val, initial, n);
316          if (!cholmod_sdmult(Ztcp, 1, one, one, chb, chv, &c))      else AZERO(val, n);
317        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"));              error(_("Error return from sdmult"));
320          Free(chv); Free(chb); cholmod_free_sparse(&Ztcp, &c);      Free(chv); Free(chb); Free(Zt);
     }  
321      return val;      return val;
322  }  }
323    
# Line 844  Line 743 
743   */   */
744  static void  static void
745  internal_glmer_reweight(GlmerStruct GS) {  internal_glmer_reweight(GlmerStruct GS) {
746      SEXP dmu_deta, var,      SEXP dmu_deta, var;
747          mZt = GET_SLOT(GS->mer, Matrix_ZtSym);      int i;
748      int *Zp = INTEGER(GET_SLOT(mZt, Matrix_pSym));      double *w = REAL(GET_SLOT(GS->mer, Matrix_wtsSym)),
749      int i, j,          *y = REAL(GET_SLOT(GS->mer, Matrix_ySym)),
750          n = LENGTH(GET_SLOT(GS->mer, Matrix_ySym)),          *z = REAL(GET_SLOT(GS->mer, Matrix_wrkresSym));
         p = LENGTH(GET_SLOT(GS->mer, Matrix_rXySym));  
     double *mX = REAL(GET_SLOT(GS->mer, Matrix_XSym)),  
         *my = REAL(GET_SLOT(GS->mer, Matrix_ySym)),  
         *mZx = REAL(GET_SLOT(mZt, Matrix_xSym));  
751    
752                                  /* reweight mer */                                  /* reweight mer */
753      eval_check_store(GS->linkinv, GS->rho, GS->mu);      eval_check_store(GS->linkinv, GS->rho, GS->mu);
# Line 861  Line 756 
756      var = PROTECT(eval_check(GS->var, GS->rho,      var = PROTECT(eval_check(GS->var, GS->rho,
757                               REALSXP, GS->n));                               REALSXP, GS->n));
758      for (i = 0; i < GS->n; i++) {      for (i = 0; i < GS->n; i++) {
759          GS->w[i] = GS->wts[i] *          w[i] = GS->wts[i] *
760              REAL(dmu_deta)[i]/sqrt(REAL(var)[i]);              REAL(dmu_deta)[i]/sqrt(REAL(var)[i]);
761          GS->z[i] = REAL(GS->eta)[i] - GS->offset[i] +          z[i] = REAL(GS->eta)[i] - GS->offset[i] +
762              (GS->y[i] - REAL(GS->mu)[i])/REAL(dmu_deta)[i];              (y[i] - REAL(GS->mu)[i])/REAL(dmu_deta)[i];
763      }      }
764      UNPROTECT(2);      UNPROTECT(2);
     for (i = 0; i < n; i++) my[i] = GS->z[i] * GS->w[i];  
     for (j = 0; j < p; j++)  
         for (i = 0; i < n; i++) mX[i + j * n] =  
                                     GS->X[i + j * n] * GS->w[i];  
     for (j = 0; j < n; j++) {  
         int i2 = Zp[j + 1];  
         for (i = Zp[j]; i < i2; i++) mZx[i] = GS->Zt[i] * GS->w[j];  
     }  
765      mer_update_ZXy(GS->mer);      mer_update_ZXy(GS->mer);
766  }  }
767    
# Line 905  Line 792 
792  }  }
793    
794  /**  /**
  * Evaluate the quadratic form (x-mn)'A'A(x-mn) from the multivariate  
  * normal kernel.  
  *  
  * @param n dimension of random variate  
  * @param mn mean vector  
  * @param a upper Cholesky factor of the precision matrix  
  * @param lda leading dimension of A  
  * @param x vector at which to evaluate the kernel  
  *  
  * @return value of the normal kernel  
  */  
 static double  
 normal_kernel(int n, const double mn[],  
               const double a[], int lda, const double x[])  
 {  
     int i, ione = 1;  
     double *tmp = Calloc(n, double), ans;  
   
     for (i = 0; i < n; i++) tmp[i] = x[i] - mn[i];  
     F77_CALL(dtrmv)("U", "N", "N", &n, a, &lda, tmp, &ione);  
     for (i = 0, ans = 0; i < n; i++) ans += tmp[i] * tmp[i];  
     Free(tmp);  
     return ans;  
 }  
   
 /**  
  * Evaluate the conditional deviance for a given set of fixed effects.  
  *  
  * @param GS Pointer to a GlmerStruct  
  * @param fixed value of the fixed effects  
  *  
  * @return conditional deviance  
  */  
 static double  
 fixed_effects_deviance(GlmerStruct GS, const double fixed[])  
 {  
     SEXP devs;  
     int i, ione = 1;  
     double ans, one = 1, zero = 0;  
   
     F77_CALL(dgemv)("N", &(GS->n), &(GS->p), &one,  
                     GS->X, &(GS->n), fixed,  
                     &ione, &zero, REAL(GS->eta), &ione);  
                                 /* add in random effects and offset */  
     vecIncrement(REAL(GS->eta), GS->off, GS->n);  
     eval_check_store(GS->linkinv, GS->rho, GS->mu);  
     devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));  
     for (i = 0, ans = 0; i < GS->n; i++) ans += REAL(devs)[i];  
     UNPROTECT(1);  
     return ans;  
 }  
   
 /**  
795   * Update the ranef slot assuming that the fixef, rZy, RZX and L slots   * Update the ranef slot assuming that the fixef, rZy, RZX and L slots
796   * are up to date.   * are up to date.
797   *   *
# Line 1027  Line 861 
861  }  }
862    
863  /**  /**
864   * Establish off, the effective offset for the fixed effects, and   * Iterate to determine the conditional modes of the random effects.
  * iterate to determine the conditional modes of the random effects.  
865   *   *
866   * @param GS a GlmerStruct object   * @param GS a GlmerStruct object
867   * @param fixed vector of fixed effects   * @param fixed vector of fixed effects
# Line 1049  Line 882 
882      Memcpy(REAL(fixef), fixed, LENGTH(fixef));      Memcpy(REAL(fixef), fixed, LENGTH(fixef));
883      internal_mer_Zfactor(GS->mer, L);      internal_mer_Zfactor(GS->mer, L);
884      internal_mer_ranef(GS->mer);      internal_mer_ranef(GS->mer);
885      internal_mer_fitted(GS->mer, GS->X, GS->Zt, GS->offset, REAL(GS->eta));      internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
886      Memcpy(GS->etaold, REAL(GS->eta), GS->n);      Memcpy(GS->etaold, REAL(GS->eta), GS->n);
887    
888      for (i = 0; i < GS->maxiter && crit > GS->tol; i++) {      for (i = 0; i < GS->maxiter && crit > GS->tol; i++) {
889          internal_glmer_reweight(GS);          internal_glmer_reweight(GS);
890          internal_mer_Zfactor(GS->mer, L);          internal_mer_Zfactor(GS->mer, L);
891          internal_mer_ranef(GS->mer);          internal_mer_ranef(GS->mer);
892          internal_mer_fitted(GS->mer, GS->X, GS->Zt, GS->offset,          internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
                             REAL(GS->eta));  
893          crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);          crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);
894      }      }
895      Free(L);      Free(L);
# Line 1170  Line 1002 
1002      }      }
1003  }  }
1004    
 /* FIXME: Separate the calculation of the offset random effects from  
  * the calculation of the deviance contributions. */  
   
 /**  
  * Determine the deviance components associated with each of the  
  * levels of a grouping factor at the conditional modes or a value  
  * offset from the conditional modes by delb.  
  *  
  * @param GS pointer to a GlmerStruct  
  * @param b conditional modes of the random effects  
  * @param Gp group pointers  
  * @param nc number of columns in the model matrix for the kth  
  * grouping factor  
  * @param k index (0-based) of the grouping factor  
  * @param delb vector of length nc giving the changes in the  
  * orthonormalized random effects  
  * @param OmgFac Cholesky factor of the inverse of the penalty matrix  
  * for this grouping factor  
  * @param bVfac 3-dimensional array holding the factors of the  
  * conditional variance-covariance matrix of the random effects  
 FIXME: This is wrong.  It is bVar[[i]] not bVfac that is being passed.  
 This only affects the AGQ method.  
  * @param devcmp array to hold the deviance components  
  *  
  * @return devcmp  
  */  
 static double*  
 rel_dev_1(GlmerStruct GS, const double b[], int nlev, int nc, int k,  
           const double delb[], const double OmgFac[],  
           const double bVfac[], double devcmp[])  
 {  
     SEXP devs;  
     int *fv = INTEGER(VECTOR_ELT(GET_SLOT(GS->mer, Matrix_flistSym), k)),  
         i, j;  
 /*     double *bcp = (double *) NULL; */  
   
     AZERO(devcmp, nlev);  
 #if 0  
     if (delb) {  
         int ione = 1, ntot = nlev * nc;  
         double sumsq = 0;  
                                 /* copy the contents of b */  
         bcp = Memcpy(Calloc(ntot, double), b, ntot);  
         if (nc == 1) {  
             sumsq = delb[0] * delb[0];  
             for (i = 0; i < nlev; i++) b[i] += delb[0] * bVfac[i];  
         } else {  
             int ncsq = nc * nc;  
             double *tmp = Calloc(nc, double);  
             for (i = 0; i < nlev; i++) {  
                 Memcpy(tmp, delb, nc);  
                 F77_CALL(dtrmv)("U", "N", "N", &nc, &(bVfac[i * ncsq]),  
                                 &nc, tmp, &ione);  
                 for (j = 0; j < nc; j++) b[i + j * nc] = tmp[j];  
             }  
                                 /* sum of squares of delb */  
             for (j = 0; j < nc; j++) sumsq += delb[j] * delb[j];  
         }  
         for (i = 0; i < nlev; i++) devcmp[i] = -sumsq;  
     }  
 #endif  
     internal_mer_fitted(  
         GS->mer, (double *) NULL,  
         REAL(GET_SLOT(GET_SLOT(GS->mer, Matrix_ZtSym), Matrix_xSym)),  
         GS->off, REAL(GS->eta));  
     eval_check_store(GS->linkinv, GS->rho, GS->mu);  
     devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));  
     for (i = 0; i < GS->n; i++)  
         devcmp[fv[i] - 1] += REAL(devs)[i];  
     UNPROTECT(1);  
     if (nc == 1) {  
         for (i = 0; i < nlev; i++) {  
             double tmp = *OmgFac * b[i];  
             devcmp[i] += tmp * tmp;  
         }  
     } else {  
         double *tmp = Calloc(nc, double);  
         int ione = 1;  
   
         for (i = 0; i < nlev; i++) {  
             for (j = 0; j < nc; j++) tmp[j] = b[i + j * nlev];  
             F77_CALL(dtrmv)("U", "N", "N", &nc, OmgFac, &nc,  
                             tmp, &ione);  
             for (j = 0; j < nc; j++)  
                 devcmp[i] += tmp[j] * tmp[j];  
         }  
     }  
 #if 0  
     if (delb) {  
         Memcpy(b, bcp, ntot);  
         Free(bcp);  
     }  
 #endif  
     return devcmp;  
 }  
   
   
1005  /**  /**
1006   * Evaluate the conditional deviance for the stored random effects.   * Evaluate the conditional deviance for the stored random effects.
1007   *   *
# Line 1282  Line 1017 
1017      int i;      int i;
1018      double ans;      double ans;
1019    
1020      internal_mer_fitted(GS->mer, GS->X, GS->Zt, GS->offset, REAL(GS->eta));      internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
1021      eval_check_store(GS->linkinv, GS->rho, GS->mu);      eval_check_store(GS->linkinv, GS->rho, GS->mu);
1022      devs = PROTECT(eval_check(GS->dev_resids, GS->rho, REALSXP, GS->n));      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];      for (i = 0, ans = 0; i < GS->n; i++) ans += REAL(devs)[i];
# Line 1290  Line 1025 
1025      return ans;      return ans;
1026  }  }
1027    
1028  static void  /* Externally accessible functions */
 internal_glmer_ranef_update(GlmerStruct GS, SEXP b)  
 {  
     SEXP bhat, bprop = PROTECT(duplicate(b)),  
         bVar = GET_SLOT(GS->mer, Matrix_bVarSym),  
         Omega = GET_SLOT(GS->mer, Matrix_OmegaSym);  
     int i, ione = 1, j, k;  
 /*     double *b = REAL(GET_SLOT(GS->mer, Matrix_ranefSym); */  
     double devr, one = 1;  
   
     bhat = R_NilValue;  
                                 /* subtract deviance at b */  
 /*     devr = -random_effects_deviance(GS, b); */  
      for (i = 0; i < GS->nf; i++) {  
         SEXP Bi = VECTOR_ELT(b, i);  
         int *dims = INTEGER(getAttrib(Bi, R_DimSymbol));  
         int nlev = dims[0], nci = dims[1];  
         int ncisqr = nci * nci, ntot = nlev * nci;  
         double *bcopy = Calloc(ntot, double),  
             *bi = REAL(Bi),  
             *bhati = REAL(VECTOR_ELT(bhat, i)),  
             *bpropi = REAL(VECTOR_ELT(bprop, i)),  
             *bVari = REAL(VECTOR_ELT(bVar, i)),  
             *chol = Calloc(ncisqr, double),  
             *delta = Calloc(nci, double),  
             *omgfac = Memcpy(Calloc(ncisqr, double),  
                              REAL(VECTOR_ELT(Omega, i)),  
                              ncisqr);  
   
                                 /* subtract quadratic form in */  
                                 /* Omega[[i]] at b  */  
         F77_CALL(dpotrf)("U", &nci, omgfac, &nci, &j);  
         if (j)  
             error(_("Leading %d minor of Omega[[%d]] not positive definite"),  
                        j, i + 1);  
         Memcpy(bcopy, bi, ntot);  
         F77_CALL(dtrmm)("R", "U", "T", "N", &nlev, &nci, &one,  
                         omgfac, &nci, bcopy, &nlev);  
         for (k = 0; k < ntot; k++) devr -= bcopy[k] * bcopy[k];  
                                 /* form bprop and proposal density */  
         for (k = 0; k < nlev; k++) {  
                                 /* proposal density at b */  
             for (j = 0; j < nci; j++)  
                 delta[j] = bi[k + j * nlev] - bhati[k + j * nlev];  
             Memcpy(chol, &(bVari[k * ncisqr]), ncisqr);  
             F77_CALL(dpotrf)("U", &nci, chol, &nci, &j);  
             if (j)  
                 error(_("Leading %d minor of bVar[[%d]][,,%d] not positive definite"),  
                        j, i + 1, k + 1);  
             F77_CALL(dtrsv)("U", "T", "N", &nci, chol, &nci,  
                             delta, &ione);  
             for (j = 0; j < nci; j++) {  
                 double nrm = norm_rand(); /* proposal deviate */  
                 devr += delta[j] * delta[j] - nrm * nrm;  
                 delta[j] = nrm;  
             }  
                                 /* scale by Cholesky inverse */  
             F77_CALL(dtrmv)("U", "T", "N", &nci, chol, &nci,  
                             delta, &ione);  
                                 /* add mean */  
             for (j = 0; j < nci; j++)  
                 bpropi[k + j * nlev] = bhati[k + j * nlev] + delta[j];  
         }  
                                 /* add quadratic form in */  
                                 /* Omega[[i]] at bprop  */  
         Memcpy(bcopy, bpropi, ntot);  
         F77_CALL(dtrmm)("R", "U", "T", "N", &nlev, &nci, &one,  
                         omgfac, &nci, bcopy, &nlev);  
         for (k = 0; k < ntot; k++) devr += bcopy[k] * bcopy[k];  
   
         Free(bcopy); Free(chol); Free(delta); Free(omgfac);  
      }  
                                 /* add deviance at bprop */  
 /*     devr += random_effects_deviance(GS, bprop); */  
   
     if (unif_rand() < exp(-0.5 * devr))  
         for (i = 0; i < GS->nf; i++) { /* copy each face of b */  
             SEXP Bi = VECTOR_ELT(b, i);  
             int *dims = INTEGER(getAttrib(Bi, R_DimSymbol));  
   
             Memcpy(REAL(Bi), REAL(VECTOR_ELT(bprop, i)),  
                    dims[0] * dims[1]);  
         }  
   
     if (asLogical(Matrix_getElement(GS->cv, "msVerbose"))) {  
         double *b0 = REAL(VECTOR_ELT(bprop, 0));  
         Rprintf("%5.3f:", exp(-0.5 * devr));  
         for (k = 0; k < 5; k++) Rprintf("%#10g ", b0[k]);  
         Rprintf("\n");  
     }  
   
     UNPROTECT(2);  
 }  
   
   
 /**  
  * Determine the conditional modes and the conditional variance of the  
  * fixed effects given the data and the current random effects.  
  * Create a Metropolis-Hasting proposal step from the multivariate  
  * normal density, determine the acceptance probability and, if the  
  * step is to be accepted, overwrite the contents of fixed with the  
  * new contents.  
  *  
  * @param GS a GlmerStruct  
  * @param b list of random effects  
  * @param fixed current value of the fixed effects  
  *  
  * @return updated value of the fixed effects  
  */  
 static double *  
 internal_glmer_fixef_update(GlmerStruct GS, SEXP b,  
                             double fixed[])  
 {  
     SEXP dmu_deta, var;  
     int i, ione = 1, it, j, lwork = -1;  
     double *ans = Calloc(GS->p, double), /* proposal point */  
         *md = Calloc(GS->p, double), /* conditional modes */  
         *w = Calloc(GS->n, double), *work,  
         *wtd = Calloc(GS->n * GS->p, double),  
         *z = Calloc(GS->n, double),  
         crit, devr, one = 1, tmp, zero = 0;  
   
     if (!isNewList(b) || LENGTH(b) != GS->nf)  
         error(_("%s must be a %s of length %d"), "b", "list", GS->nf);  
     for (i = 0; i < GS->nf; i++) {  
         SEXP bi = VECTOR_ELT(b, i);  
         if (!isReal(bi) || !isMatrix(bi))  
             error(_("b[[%d]] must be a numeric matrix"), i);  
     }  
     AZERO(z, GS->n);            /* -Wall */  
     Memcpy(md, fixed, GS->p);  
                                 /* calculate optimal size of work array */  
     F77_CALL(dgels)("N", &(GS->n), &(GS->p), &ione, wtd, &(GS->n),  
                     z,  &(GS->n), &tmp, &lwork, &j);  
     if (j)                      /* shouldn't happen */  
         error(_("%s returned error code %d"), "dgels", j);  
     lwork = (int) tmp;  
     work = Calloc(lwork, double);  
   
     AZERO(GS->off, GS->n); /* fitted values from random effects */  
 /*     fitted_ranef(GET_SLOT(GS->mer, Matrix_flistSym), GS->unwtd, b, */  
 /*               INTEGER(GET_SLOT(GS->mer, Matrix_ncSym)), GS->off); */  
     for (i = 0; i < GS->n; i++)  
         (GS->etaold)[i] = ((GS->off)[i] += (GS->offset)[i]);  
   
     for (it = 0, crit = GS->tol + 1;  
          it < GS->maxiter && crit > GS->tol; it++) {  
                                 /* fitted values from current beta */  
         F77_CALL(dgemv)("N", &(GS->n), &(GS->p), &one,  
                         GS->X, &(GS->n), md,  
                         &ione, &zero, REAL(GS->eta), &ione);  
                                 /* add in random effects and offset */  
         vecIncrement(REAL(GS->eta), (GS->off), GS->n);  
                                 /* check for convergence */  
         crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);  
                                 /* obtain mu, dmu_deta, var */  
         eval_check_store(GS->linkinv, GS->rho, GS->mu);  
         dmu_deta = PROTECT(eval_check(GS->mu_eta, GS->rho,  
                                       REALSXP, GS->n));  
         var = PROTECT(eval_check(GS->var, GS->rho, REALSXP, GS->n));  
                                 /* calculate weights and working residual */  
         for (i = 0; i < GS->n; i++) {  
             w[i] = GS->wts[i] * REAL(dmu_deta)[i]/sqrt(REAL(var)[i]);  
             z[i] = w[i] * (REAL(GS->eta)[i] - (GS->off)[i] +  
                            ((GS->y)[i] - REAL(GS->mu)[i]) /  
                            REAL(dmu_deta)[i]);  
         }  
         UNPROTECT(2);  
                                 /* weighted copy of the model matrix */  
         for (j = 0; j < GS->p; j++)  
             for (i = 0; i < GS->n; i++)  
                 wtd[i + j * GS->n] = GS->X[i + j * GS->n] * w[i];  
                                 /* weighted least squares solution */  
         F77_CALL(dgels)("N", &(GS->n), &(GS->p), &ione, wtd, &(GS->n),  
                         z, &(GS->n), work, &lwork, &j);  
         if (j) error(_("%s returned error code %d"), "dgels", j);  
         Memcpy(md, z, GS->p);  
     }  
                                 /* wtd contains the Cholesky factor of  
                                  * the precision matrix */  
     devr = normal_kernel(GS->p, md, wtd, GS->n, fixed);  
     devr -= fixed_effects_deviance(GS, fixed);  
     for (i = 0; i < GS->p; i++) {  
         double var = norm_rand();  
         ans[i] = var;  
         devr -= var * var;  
     }  
     F77_CALL(dtrsv)("U", "N", "N", &(GS->p), wtd, &(GS->n), ans, &ione);  
     for (i = 0; i < GS->p; i++) ans[i] += md[i];  
     devr += fixed_effects_deviance(GS, ans);  
     crit = exp(-0.5 * devr);    /* acceptance probability */  
     tmp = unif_rand();  
     if (asLogical(Matrix_getElement(GS->cv, "msVerbose"))) {  
         Rprintf("%5.3f: ", crit);  
         for (j = 0; j < GS->p; j++) Rprintf("%#10g ", ans[j]);  
         Rprintf("\n");  
     }  
     if (tmp < crit) Memcpy(fixed, ans, GS->p);  
     Free(ans); Free(md); Free(w);  
     Free(work); Free(wtd); Free(z);  
     return fixed;  
 }  
   
 /* Gauss-Hermite Quadrature x positions and weights */  
 static const double  
     GHQ_x1[1] = {0},  
     GHQ_w1[1] = {1},  
     GHQ_x2[1] = {1},  
     GHQ_w2[1] = {0.5},  
     GHQ_x3[2] = {1.7320507779261, 0},  
     GHQ_w3[2] = {0.166666666666667, 0.666666666666667},  
     GHQ_x4[2] = {2.3344141783872, 0.74196377160456},  
     GHQ_w4[2] = {0.0458758533899086, 0.454124131589555},  
     GHQ_x5[3] = {2.85696996497785, 1.35562615677371, 0},  
     GHQ_w5[3] = {0.0112574109895360, 0.222075915334214,  
                  0.533333317311434},  
     GHQ_x6[3] = {3.32425737665988, 1.88917584542184,  
                  0.61670657963811},  
     GHQ_w6[3] = {0.00255578432527774, 0.0886157433798025,  
                  0.408828457274383},  
     GHQ_x7[4] = {3.7504396535397, 2.36675937022918,  
                  1.15440537498316, 0},  
     GHQ_w7[4] = {0.000548268839501628, 0.0307571230436095,  
                  0.240123171391455, 0.457142843409801},  
     GHQ_x8[4] = {4.14454711519499, 2.80248581332504,  
                  1.63651901442728, 0.539079802125417},  
     GHQ_w8[4] = {0.000112614534992306, 0.0096352198313359,  
                  0.117239904139746, 0.373012246473389},  
     GHQ_x9[5] = {4.51274578616743, 3.20542894799789,  
                  2.07684794313409, 1.02325564627686, 0},  
     GHQ_w9[5] = {2.23458433364535e-05, 0.00278914123744297,  
                  0.0499164052656755, 0.244097495561989,  
                  0.406349194142045},  
     GHQ_x10[5] = {4.85946274516615, 3.58182342225163,  
                   2.48432579912153, 1.46598906930182,  
                   0.484935699216176},  
     GHQ_w10[5] = {4.31065250122166e-06, 0.000758070911538954,  
                   0.0191115799266379, 0.135483698910192,  
                   0.344642324578594},  
     GHQ_x11[6] = {5.18800113558601, 3.93616653976536,  
                   2.86512311160915, 1.87603498804787,  
                   0.928868981484148, 0},  
     GHQ_w11[6] = {8.12184954622583e-07, 0.000195671924393029,  
                   0.0067202850336527, 0.066138744084179,  
                   0.242240292596812, 0.36940835831095};  
   
 static const double  
     *GHQ_x[12] = {(double *) NULL, GHQ_x1, GHQ_x2, GHQ_x3, GHQ_x4,  
                   GHQ_x5, GHQ_x6, GHQ_x7, GHQ_x8, GHQ_x9, GHQ_x10,  
                   GHQ_x11},  
     *GHQ_w[12] = {(double *) NULL, GHQ_w1, GHQ_w2, GHQ_w3, GHQ_w4,  
                   GHQ_w5, GHQ_w6, GHQ_w7, GHQ_w8, GHQ_w9, GHQ_w10,  
                   GHQ_w11};  
   
 /* Externally accessible functions */  
1029    
1030  /**  /**
1031   * Simulate a sample of random matrices from a Wishart distribution   * Simulate a sample of random matrices from a Wishart distribution
# Line 1613  Line 1094 
1094          internal_ECMEsteps(GS->mer, i ? 2 : GS->niterEM,          internal_ECMEsteps(GS->mer, i ? 2 : GS->niterEM,
1095                             GS->EMverbose);                             GS->EMverbose);
1096          eval(GS->LMEopt, GS->rho);          eval(GS->LMEopt, GS->rho);
1097          internal_mer_fitted(GS->mer, GS->X, GS->Zt, GS->offset, REAL(GS->eta));          internal_mer_fitted(GS->mer, GS->offset, REAL(GS->eta));
1098          crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);          crit = conv_crit(GS->etaold, REAL(GS->eta), GS->n);
1099      }      }
1100      if (crit > GS->tol)      if (crit > GS->tol)
# Line 1622  Line 1103 
1103      return R_NilValue;      return R_NilValue;
1104  }  }
1105    
 #if 0  
 SEXP glmer_bhat(SEXP pars, SEXP GSp)  
 {  
     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);  
   
     if (!isReal(pars) || LENGTH(pars) != GS->npar)  
         error(_("`%s' must be a numeric vector of length %d"),  
               "pars", GS->npar);  
     if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))  
         warning(_("internal_bhat did not converge"));  
     return R_NilValue;  
 }  
 #endif  
   
1106  /**  /**
1107   * Compute the Laplace approximation to the deviance.   * Compute the Laplace approximation to the deviance.
1108   *   *
# Line 1666  Line 1133 
1133  }  }
1134    
1135  /**  /**
1136   * Compute the approximation to the deviance using adaptive   * Release the storage for a GlmerStruct
  * Gauss-Hermite quadrature (AGQ).  When nAGQ == 1 this is the Laplace  
  * approximation.  
  *  
  * @param pars pointer to a numeric vector of parameters  
  * @param GSp pointer to a GlmerStruct object  
  * @param nAGQp pointer to a scalar integer representing the number of  
  * points in AGQ to use  
  *  
  * @return the approximation to the deviance as computed using AGQ  
  */  
 SEXP glmer_devAGQ(SEXP pars, SEXP GSp, SEXP nAGQp)  
 {  
     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);  
     SEXP Omega = GET_SLOT(GS->mer, Matrix_OmegaSym),  
         bVar = GET_SLOT(GS->mer, Matrix_bVarSym);  
     int i, j, k, nAGQ = asInteger(nAGQp);  
     int n2 = (nAGQ + 1)/2,  
         *Gp = INTEGER(GET_SLOT(GS->mer, Matrix_GpSym)),  
         *nc = INTEGER(GET_SLOT(GS->mer, Matrix_ncSym));  
     double *f0, LaplaceDev = 0, AGQadjst = 0,  
         *bhat = REAL(GET_SLOT(GS->mer, Matrix_ranefSym));  
   
     if (!isReal(pars) || LENGTH(pars) != GS->npar)  
         error(_("`%s' must be a numeric vector of length %d"),  
               "pars", GS->npar);  
     if (GS->nf > 1 && nAGQ > 1) {  
         warning(_("AGQ not available for multiple grouping factors - using Laplace"));  
         nAGQ = 1;  
     }  
     if (!internal_bhat(GS, REAL(pars), REAL(pars) + (GS->p)))  
         return ScalarReal(DBL_MAX);  
   
     for (i = 0; i < GS->nf; i++) {  
         int nci = nc[i];  
         int ncip1 = nci + 1, ncisqr = nci * nci,  
             nlev = (Gp[i + 1] - Gp[i])/nci;  
         double *omgf = REAL(GET_SLOT(dpoMatrix_chol(VECTOR_ELT(Omega, i)), Matrix_xSym)),  
             *bVi = Memcpy(Calloc(ncisqr * nlev, double),  
                            REAL(VECTOR_ELT(bVar, i)), ncisqr * nlev);  
   
         for (j = 0; j < nci; j++) { /* nlev * logDet(Omega_i) */  
             LaplaceDev += 2 * nlev * log(omgf[j * ncip1]);  
         }  
         for (k = 0; k < nlev; k++) {  
             double *bVik = bVi + k * ncisqr;  
             F77_CALL(dpotrf)("U", &nci, bVik, &nci, &j);  
             if (j)  
                 error(_("Leading %d minor of bVar[[%d]][,,%d] not positive definite"),  
                       j, i + 1, k + 1);  
             for (j = 0; j < nci; j++) LaplaceDev -= 2 * log(bVik[j * ncip1]);  
         }  
   
         f0 = Calloc(nlev, double);  
         rel_dev_1(GS, bhat, nlev, nci, i, (double *) NULL,  
                   omgf, bVi, f0);  
         for (k = 0; k < nlev; k++) LaplaceDev += f0[k];  
         if (nAGQ > 1) {  
             double *fx = Calloc(nlev, double),  
                 *rellik = Calloc(nlev, double),  
                 *delb = Calloc(nci, double);  
   
             if (nci > 1) error(_("code not yet written"));  
             AZERO(rellik, nlev);        /* zero accumulator */  
             for (k = 0; k < n2; k++) {  
                 delb[0] = GHQ_x[nAGQ][k];  
                 if (delb[0]) {  
                     rel_dev_1(GS, bhat, nlev, nci, i, delb,  
                               omgf, bVi, fx);  
                     for (j = 0; j < nlev; j++) {  
                         rellik[j] += GHQ_w[nAGQ][k] *  
                             exp(-(fx[j] - f0[j])/2);  
                     }  
                     delb[0] *= -1;  
                     rel_dev_1(GS, bhat, nlev, nci, i, delb,  
                               omgf, bVi, fx);  
                     for (j = 0; j < nlev; j++) {  
                         rellik[j] += GHQ_w[nAGQ][k] *  
                             exp(-(fx[j] - f0[j])/2);  
                     }  
                 } else {  
                     for (j = 0; j < nlev; j++)  
                         rellik[j] += GHQ_w[nAGQ][k];  
                 }  
             }  
             for (j = 0; j < nlev; j++)  
                 AGQadjst -= 2 * log(rellik[j]);  
             Free(fx); Free(rellik);  
         }  
         Free(f0); Free(bVi);  
     }  
     return ScalarReal(LaplaceDev + AGQadjst);  
 }  
   
 /**  
  * Release the storage for a GlmerStruct  
1137   *   *
1138   * @param GSp External pointer to a  GlmerStruct   * @param GSp External pointer to a  GlmerStruct
1139   *   *
# Line 1770  Line 1142 
1142  SEXP glmer_finalize(SEXP GSp) {  SEXP glmer_finalize(SEXP GSp) {
1143      GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);      GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);
1144    
1145      Free(GS->w); Free(GS->X); Free(GS->y); Free(GS->z);      Free(GS->offset); Free(GS->wts); Free(GS->etaold);
     Free(GS->Zt); Free(GS->off); Free(GS->offset); Free(GS->wts);  
     Free(GS->etaold);  
1146      Free(GS);      Free(GS);
1147      return R_NilValue;      return R_NilValue;
1148  }  }
1149    
1150  /**  /**
  * Determine the conditional modes and the conditional variance of the  
  * fixed effects given the data and the current random effects.  
  * Create a Metropolis-Hasting proposal step from the multivariate  
  * normal density, determine the acceptance probability and return the  
  * current value or the proposed value.  
  *  
  * @param GSp pointer to a GlmerStruct  
  * @param b list of random effects  
  * @param fixed current value of the fixed effects  
  *  
  * @return updated value of the fixed effects  
  */  
 SEXP glmer_fixed_update(SEXP GSp, SEXP b, SEXP fixed)  
 {  
     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);  
   
     if (!isReal(fixed) || LENGTH(fixed) != GS->p)  
         error(_("%s must be a %s of length %d"), "fixed",  
                 "numeric vector", GS->p);  
     GetRNGstate();  
     internal_glmer_fixef_update(GS, b, REAL(fixed));  
     PutRNGstate();  
     return fixed;  
 }  
   
 /**  
1151   * Return an external pointer object to a GlmerStruct created in   * Return an external pointer object to a GlmerStruct created in
1152   * environment rho   * environment rho
1153   *   *
# Line 1824  Line 1168 
1168      y = GET_SLOT(GS->mer, Matrix_ySym);      y = GET_SLOT(GS->mer, Matrix_ySym);
1169      GS->n = LENGTH(y);      GS->n = LENGTH(y);
1170      GS->p = LENGTH(GET_SLOT(GS->mer, Matrix_rXySym));      GS->p = LENGTH(GET_SLOT(GS->mer, Matrix_rXySym));
     GS->X = Memcpy(Calloc(GS->n * GS->p, double),  
                    REAL(GET_SLOT(GS->mer, Matrix_XSym)), GS->n * GS->p);  
1171      GS->y = Memcpy(Calloc(GS->n, double), REAL(y), GS->n);      GS->y = Memcpy(Calloc(GS->n, double), REAL(y), GS->n);
     GS->w = Calloc(GS->n, double);  
     GS->z = Calloc(GS->n, double);  
1172      Ztx = GET_SLOT(GET_SLOT(GS->mer, Matrix_ZtSym), Matrix_xSym);      Ztx = GET_SLOT(GET_SLOT(GS->mer, Matrix_ZtSym), Matrix_xSym);
     GS->Zt = Memcpy(Calloc(LENGTH(Ztx), double), REAL(Ztx), LENGTH(Ztx));  
1173      GS->mu = find_and_check(rho, install("mu"), REALSXP, GS->n);      GS->mu = find_and_check(rho, install("mu"), REALSXP, GS->n);
1174      tmp = find_and_check(rho, install("offset"), REALSXP, GS->n);      tmp = find_and_check(rho, install("offset"), REALSXP, GS->n);
1175      GS->offset = Memcpy(Calloc(GS->n, double), REAL(tmp), GS->n);      GS->offset = Memcpy(Calloc(GS->n, double), REAL(tmp), GS->n);
1176      tmp = find_and_check(rho, install("wts"), REALSXP, GS->n);      tmp = find_and_check(rho, install("wts"), REALSXP, GS->n);
1177      GS->wts = Memcpy(Calloc(GS->n, double), REAL(tmp), GS->n);      GS->wts = Memcpy(Calloc(GS->n, double), REAL(tmp), GS->n);
     GS->off = Calloc(GS->n, double);  
1178      GS->etaold = Calloc(GS->n, double);      GS->etaold = Calloc(GS->n, double);
1179      GS->cv = find_and_check(rho, install("cv"), VECSXP, 0);      GS->cv = find_and_check(rho, install("cv"), VECSXP, 0);
1180      GS->niterEM = asInteger(Matrix_getElement(GS->cv, "niterEM"));      GS->niterEM = asInteger(Matrix_getElement(GS->cv, "niterEM"));
# Line 1862  Line 1200 
1200  }  }
1201    
1202  /**  /**
  * Determine the conditional modes and the conditional variance of the  
  * random effects given the data and the current fixed effects and  
  * variance components.  
  *  
  * Create a Metropolis-Hasting proposal step from a multivariate  
  * normal density centered at bhat with variance-covariance matrix  
  * from bVar, determine the acceptance probability and return the  
  * current value or the proposed value.  
  *  
  * @param GSp pointer to a GlmerStruct  
  * @param fixed pointer to a numeric vector of the fixed effects  
  * @param varc pointer to a numeric vector of the variance components  
  * @param varc pointer to current values of b  
  *  
  * @return updated b from the Metropolis-Hastings step  
  */  
 SEXP glmer_ranef_update(SEXP GSp, SEXP fixed, SEXP varc, SEXP b)  
 {  
     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSp);  
     int nvarc = GS->npar - GS->p;  
   
     if (!isReal(fixed) || LENGTH(fixed) != GS->p)  
         error(_("`%s' must be a numeric vector of length %d"),  
               "fixed", GS->p);  
     if (INTEGER(GET_SLOT(GS->mer, Matrix_ncSym))[GS->nf] > 0)  
         error(_("the mer object must be set to skip fixed effects"));  
     if (!isReal(varc) || LENGTH(varc) != nvarc)  
         error(_("`%s' must be a numeric vector of length %d"),  
               "varc", nvarc);  
   
     GetRNGstate();  
     /* Don't check for convergence failure after internal_bhat.  
      * It is determining the mean of the proposal density and  
      * does not need exact convergence. */  
     internal_bhat(GS, REAL(fixed), REAL(varc));  
     internal_glmer_ranef_update(GS, b);  
     PutRNGstate();  
   
     UNPROTECT(1);  
     return b;  
 }  
   
 /**  
1203   * Perform ECME steps for the REML or ML criterion.   * Perform ECME steps for the REML or ML criterion.
1204   *   *
1205   * @param x pointer to an mer object   * @param x pointer to an mer object
# Line 2112  Line 1407 
1407  }  }
1408    
1409  /**  /**
  * Create a Markov Chain Monte Carlo sample from a fitted generalized  
  * linear mixed model  
  *  
  * @param GSpt External pointer to a GlmerStruct  
  * @param b Conditional modes of the random effects at the parameter  
  * estimates  
  * @param fixed Estimates of the fixed effects  
  * @param varc Estimates of the variance components  
  * @param savebp Logical indicator of whether or not to save the  
  * random effects in the MCMC sample  
  * @param nsampp Integer value of the number of samples to generate  
  *  
  * @return  
  */  
 SEXP  
 glmer_MCMCsamp(SEXP GSpt, SEXP b, SEXP fixed, SEXP varc,  
                SEXP savebp, SEXP nsampp)  
 {  
     GlmerStruct GS = (GlmerStruct) R_ExternalPtrAddr(GSpt);  
     int i, j, nf = LENGTH(b), nsamp = asInteger(nsampp),  
         p = LENGTH(fixed), q = LENGTH(varc),  
         saveb = asLogical(savebp);  
     int *mcpar, nc = p + q;  
     SEXP ans, mcparSym = install("mcpar");  
   
     if (nsamp <= 0) nsamp = 1;  
     nc = p + q;  
     if (saveb)  
         for (i = 0; i < nf; i++) {  
             int *dims = INTEGER(getAttrib(VECTOR_ELT(b, i),  
                                           R_DimSymbol));  
             nc += dims[0] * dims[1];  
         }  
     ans = PROTECT(allocMatrix(REALSXP, nsamp, nc));  
     GetRNGstate();  
     for (i = 0; i < nsamp; i++) {  
         internal_glmer_fixef_update(GS, b, REAL(fixed));  
         internal_bhat(GS, REAL(fixed), REAL(varc));  
         internal_glmer_ranef_update(GS, b);  
 /*      internal_Omega_update(GS->mer, b); */  
         internal_mer_coef(GS->mer, 2, REAL(varc));  
         for (j = 0; j < p; j++)  
             REAL(ans)[i + j * nsamp] = REAL(fixed)[j];  
         for (j = 0; j < q; j++)  
             REAL(ans)[i + (j + p) * nsamp] = REAL(varc)[j];  
         if (saveb) {  
             int base = p + q, k;  
             for (k = 0; k < nf; k++) {  
                 SEXP bk = VECTOR_ELT(b, k);  
                 int *dims = INTEGER(getAttrib(bk, R_DimSymbol));  
                 int klen = dims[0] * dims[1];  
   
                 for (j = 0; j < klen; j++)  
                     REAL(ans)[i + (j + base) * nsamp] = REAL(bk)[j];  
                 base += klen;  
             }  
         }  
     }  
     PutRNGstate();  
     UNPROTECT(1);  
                                 /* set (S3) class and mcpar attribute */  
     setAttrib(ans, R_ClassSymbol, mkString("mcmc"));  
     setAttrib(ans, mcparSym, allocVector(INTSXP, 3));  
     mcpar = INTEGER(getAttrib(ans, mcparSym));  
     mcpar[0] = mcpar[2] = 1;  
     mcpar[1] = nsamp;  
   
     return ans;  
 }  
   
 /**  
1410   * Create an mer object from a list of grouping factors and a list of model   * Create an mer object from a list of grouping factors and a list of model
1411   * matrices.   * matrices.
1412   *   *
# Line 2381  Line 1605 
1605      return x;      return x;
1606  }  }
1607    
   
1608  /**  /**
1609   * Return L as a dtCMatrix object   * Return L as a dtCMatrix object
1610   *   *
# Line 2524  Line 1747 
1747   * @return pointer to a numeric array of fitted values   * @return pointer to a numeric array of fitted values
1748   */   */
1749    
1750  SEXP mer_fitted(SEXP x, SEXP useFe, SEXP useRe)  SEXP mer_fitted(SEXP x)
1751  {  {
1752      int n = LENGTH(GET_SLOT(x, Matrix_ySym));      int n = LENGTH(GET_SLOT(x, Matrix_ySym));
1753      SEXP ans = PROTECT(allocVector(REALSXP, n));      SEXP ans = PROTECT(allocVector(REALSXP, n));
1754    
1755      internal_mer_fitted(x,      internal_mer_fitted(x, (double*) NULL, REAL(ans));
                         (asLogical(useFe)  
                          ? REAL(GET_SLOT(x, Matrix_XSym))  
                          : (double *) NULL),  
                         (asLogical(useRe)  
                          ? REAL(GET_SLOT(GET_SLOT(x, Matrix_ZtSym), Matrix_xSym))  
                          : (double *) NULL),  
                         (double *) NULL, REAL(ans));  
1756      UNPROTECT(1);      UNPROTECT(1);
1757      return ans;      return ans;
1758  }  }
# Line 2789  Line 2005 
2005  }  }
2006    
2007  /**  /**
  * Return the permutation of the columns of Z as a pMatrix object  
  *  
  * @param x pointer to an mer object  
  *  
  * @return the permutation as an pMatrix object  
  */  
 SEXP mer_pMatrix(SEXP x)  
 {  
     cholmod_factor *L = as_cholmod_factor(GET_SLOT(x, Matrix_LSym));  
     SEXP ans = PROTECT(NEW_OBJECT(MAKE_CLASS("pMatrix")));  
     int *dims = INTEGER(ALLOC_SLOT(ans, Matrix_DimSym, INTSXP, 2)),  
         *perm = INTEGER(ALLOC_SLOT(ans, Matrix_permSym, INTSXP, L->n)), i;  
   
     dims[0] = dims[1] = (int) L->n;  
     for (i = 0; i < (int) L->n; i++) perm[i] = ((int *)(L->Perm))[i] + 1;  
     Free(L);  
     UNPROTECT(1);  
     return ans;  
 }  
   
 /**  
2008   * Extract the conditional modes of the random effects.   * Extract the conditional modes of the random effects.
2009   *   *
2010   * @param x Pointer to an mer object   * @param x Pointer to an mer object
# Line 3086  Line 2281 
2281      UNPROTECT(1);      UNPROTECT(1);
2282      return ans;      return ans;
2283  }  }
2284    
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

Legend:
Removed from v.1161  
changed lines
  Added in v.1162

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