/* Procedures and function for generalized linear mixed and linear
   mixed models */

fn DIGAM(b) = ln(b) - 0.5./b - (1/12)./(b.*b) + (1./120).*b^(-4)
 - (1/252).*b^(-6) + (1/240).*b^(-8) - (1/132)*b^(-10) + (691/2730/12)*b^(-12);

fn TRIGAM(b) = 1./b + 1./(2*b.*b) + (1./6)*b^(-3)
 - (1./30)*b^(-5) + (1./42)*b^(-7) - (1/30)*b^(-9) + (5/66)*b^(-11);

PROC(4) = ML(s,r,z,q,x,y,test,lflag);
  cls;
  output on;"********** EM algorithm for Linear Mixed Model *************";
  output off;

  if r /= rows(s);
    "s0 is wrong length";
  endif;

  /* setup indexes for random effects */
  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  /* initialize */
  n = rows(y);
  b = zeros(qs,1);
  s0 = vcx(y);
  a = invpd(x'x)'x'y;
  p = cols(x);
  nit = 0;
  iter = {};
  locate 3,1;"iter    llk     s0     a ... s  ";

  /* start of EM cycles.  EM cycles until converge = 1 */
  converge = 0;
  do until converge;

    olpv = pv;

    /* build matrix D and D inverse */
    I = eye(qs);
    d = {};
    j = 1;
    do until j == r + 1;
      t = seqa(qidx[j,1],1,q[j]);
      d = d | s[j]*ones(q[j],1);
      j = j + 1;
    endo;
    din = diagrv(I,1./d);

    /* compute conditional mean (m) */
    m = X*a + Z*b;

    /* Compute b(p) */

    vin = invpd(z'*z/s0 + din);
    b = vin*Z'(y - X*a)/s0;

    /* M-step for alp */
    a = invpd(x'x)*x'(y - Z*b);
    s0 = (1/n)*( (y-m)'(y-m) + tr(vin*Z'Z) );

    /* M-step for s */
    j = 1;
    do until j == r + 1;
      t = seqa(qidx[j,1],1,q[j]);
      s[j] = (b[t]'b[t] + tr(vin[t,t]))/q[j];
      j = j + 1;
    endo;

    /* pv is parameter vector */
    pv = a|s0|s;
    iter = iter|pv';  @ save iterates @
    converge = rows(pv) == sumc(abs(pv-olpv) .< test);

    locate 4,1;
    format 4,0;nit;;format 7,3;
    @ compute log-likehood @
    if lflag*(nit > 1);
        Hin = (1/s0)*(eye(n) - Z*vin*Z'/s0);
        t = -(1/2)*((y-X*a)'Hin*(y-X*a) - ln(det(Hin)));
    else;
      t = 0;
    endif;
    st = st|t;  @ vector of adjusted logL's@
    t;;s0;;a';;s';
    nit = nit + 1;
  endo;
  output on;
  ?;"EM MLEs for linear mixed model ... ";
  "Fixed effects estimates";format 8,4;a';
  "Residual variance";s0;
  "Variance component estimates";format 8,4;s';
  output off;

RETP(a,s,s0,b);
ENDP;

PROC(4) = REML(s,r,z,b,q,x,y,test);
  cls;output on;"****************** EM algorithm for REML *****************";
  output off;

  if r /= rows(s);
    "s0 is wrong length";
  endif;

  /* setup indexes for random effects */
  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  /* initialize */
  n = rows(y);
  b = zeros(qs,1);
  s0 = vcx(y);
  a = invpd(x'x)'x'y;
  p = cols(x);
  pv = s0|s;
  nit = 0;
  iter = {};
  ep = eye(n) - x*invpd(x'x)*x';
  locate 3,1;"iter    llk     s0     s  ";;

  /* start of EM.  EM cycles until converge = 1 */
  converge = 0;
  do until converge;

    olpv = pv;

    /* build matrix D and D inverse */
    I = eye(qs);
    d = {};
    j = 1;
    do until j == r + 1;
      t = seqa(qidx[j,1],1,q[j]);
      d = d | s[j]*ones(q[j],1);
      j = j + 1;
    endo;
    din = diagrv(I,1./d);
    D = diagrv(I,d);

    vin = invpd(z'*z/s0 + din);
    Hin = (1/s0)*(eye(n) - Z*vin*Z'/s0);
    Pm = Hin - Hin*X*invpd(X'Hin*X)*X'Hin;
    b = D*Z'Pm*y;

    /* Compute s0(p) */
    V = D - D*Z'Pm*Z*D;
    s0 = ( (y-Z*b)'ep*(y-Z*b) + tr(Z*V*Z'ep) )/(n-p);

    /* M-step for s */
    j = 1;
    do until j == r + 1;
      t = seqa(qidx[j,1],1,q[j]);
      s[j] = (b[t]'b[t] + tr(v[t,t]))/q[j];
      j = j + 1;
    endo;

    /* pv is parameter vector */
    pv = s0|s;
    iter = iter|pv';  @ save iterates @
    converge = rows(pv) == sumc(abs(pv-olpv) .< test);
    va = eigh(Pm);
    t = -(1/2)*( y'Pm*y - sumc(ln(selif(va,va.>0.001))) );
    st = st|t;  @ vector of adjusted logL's@

    locate 4,1;
    format 4,0;nit;;format 7,3;
    t;;s0;;s';
    nit = nit + 1;

  endo;
  /* M-step for alp */
  a = invpd(x'x)*x'(y - Z*b);
  output on;
  ?;"EM REML estimates for linear mixed model ... ";
  "Fixed effects estimates";format 8,4;a';
  "Residual variance ";s0;
  "Variance component estimates";format 8,4;s';
  output off;

RETP(a,s0,s,b);
ENDP;

PROC(2) = IOML(a,s,s0,b,x,z,q,names,dname);
  @ calculates information matrix for linear mixed models @
  output on;
  ?;"SE Computation for Linear Mixed Model";
  $dname;
  format 2,0;"   N = ";;rows(y);;
  p = cols(x);
  r = rows(qidx);
  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  m = X*a + Z*b;

  /* build matrix D and D inverse */
  I = eye(qs);
  d = {};
  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    d = d | s[j]*ones(q[j],1);
    j = j + 1;
  endo;
  din = diagrv(I,1./d);
  vin = invpd(z'*z/s0 + din);
  hin = eye(n)/s0 - Z*vin*Z'/(s0^2);

  ion = zeros(p+r+1,p+r+1);
  ion[1:p,1:p] = X'Hin*X;
  ion[p+1,p+1] = (1/2)*tr(Hin*Hin);

  k = 1;
  do until k == r + 1;
    t = seqa(qidx[k,1],1,q[k]);
    Ik = zeros(qs,qs);
    Ik[t,t] = eye(q[k]);

    ion[p+k+1,p+k+1] = (1/2)*tr(Z'Hin*Z*Ik*Z'Hin*Z*Ik);

    ion[p+1,p+1+k] = (1/2)*tr(Hin*Z*Ik*Z'Hin);
    ion[p+1+k,p+1] = ion[p+1,p+1+k];
    i = k + 1;
    do until i == r + 1;
      Ij = zeros(qs,qs);
      t = seqa(qidx[i,1],1,q[i]);
      Ij[t,t] = eye(q[i]);
      ion[p+1+i,p+1+k] = (1/2)*tr(Z'Hin*Z*Ik*Z'Hin*Z*Ij);

      ion[p+1+k,p+1+i] = ion[p+1+i,p+1+k];
      i = i + 1;
    endo;
    k = k + 1;
  endo;
  oldval = trapchk(1);
  trap 1,1;
  ioin = invpd(ion);
  trap oldval,1;
  if scalerr(ioin);
    "Info matrix singular";
    t = seqa(1,1,p);
    ioin = invpd(ion[t,t]);
    se = miss(ones(p+r,1),ones(p+r,1));
    se[1:p] = sqrt(diag(ioin));
  else;
    se = sqrt(diag(invpd(ion)));
  endif;

  format 8,3;
  fmt = {"*.*s" 8 8};
  t = {"*.*lf" 10 4};
  j = 1;
  do until j == 4;
    fmt = fmt|t;
    j = j + 1;
  endo;
  r = rows(s);
  names = name |"s0";
  t = se;
  t[p+1] = se[p+1]/sqrt(s0)/2;
  j = 1;
  do until j == r+ 1;
    t[p+j+1] = se[p+j+1]/sqrt(s[j])/2;
    names = names | "s"$+ftocv(j,1,0);
    j = j + 1;
  endo;

  format 8,4;

  /* transform from variance components to standard deviations */
  mask = zeros(rows(names),1)~ones(rows(names),rows(fmt)-1);
  ?;"Variable  Coefficient   SE      Est/SE ";
  ret = printfm(names~(a|sqrt(s0|s))~t~((a|sqrt(s0|s))./t),mask,fmt);?;
 "  Table reports standard deviations instead of variance component ests";

  ?;"Variance component estimates";
  "       Est        SE";
  ret = printfm(s~se[p+1:p+r],mask[p+1:p+r,2 3],fmt[2 3,.]);?;
  format 1,0;"  N iterations = ";;rows(iter);format 8,2;
  /* calc. likelihood */
  t = -(1/2)*((y-X*a)'Hin*(y-X*a) - ln(det(Hin)));
  format 8,3;"Log-likelihood = ";;t;

  output off;

RETP(ion,se);
ENDP;


PROC(2) = IOREML(a,s0,s,b,x,z,q,names,dname);
  @ calculates information matrix for REML @
  output on;
  ?;"SE Computation for REML estimation of the Linear Mixed Model";
  $dname;
  format 2,0;"   N = ";;rows(y);;
  p = cols(x);
  r = rows(qidx);
  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  m = X*a + Z*b;

  /* build matrix D and D inverse */
  I = eye(qs);
  d = {};
  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    d = d | s[j]*ones(q[j],1);
    j = j + 1;
  endo;
  din = diagrv(I,1./d);

  ion = zeros(p+r+1,p+r+1);
  vin = invpd(z'*z/s0 + din);
  Hin = (1/s0)*(eye(n) - Z*vin*Z'/s0);
  Pm = Hin - Hin*X*invpd(X'Hin*X)*X'Hin;
  ion[1:p,1:p] = X'Hin*X;
  ion[p+1,p+1] = y'Pm*Pm*Pm*y - (1/2)*tr(Pm*Pm);
  ion[p+1,p+1] = (1/2)*tr(Pm*Pm);

  k = 1;
  do until k == r + 1;

    t = seqa(qidx[k,1],1,q[k]);
    e = Z[.,t]*Z[.,t]';
    ion[p+k+1,p+k+1] = (1/2)*tr(Pm*E*Pm*E);

    ion[p+1,p+1+k] = (1/2)*tr(Pm*Pm*E);
    ion[p+1+k,p+1] = ion[p+1,p+1+k];
    i = k + 1;
    do until i == r + 1;
      u = seqa(qidx[i,1],1,q[i]);
      f = Z[.,u]*Z[.,u]';
      ion[p+k+i,p+k+k] = (1/2)*tr(Pm*E*Pm*F);

      ion[p+1+k,p+1+i] = ion[p+1+i,p+1+k];
      i = i + 1;
    endo;
    k = k + 1;
  endo;
  oldval = trapchk(1);
  trap 1,1;
  ioin = invpd(ion);
  trap oldval,1;
  if scalerr(ioin);
    "Info matrix singular";
    t = seqa(1,1,p);
    ioin = invpd(ion[t,t]);
    se = miss(ones(p+r,1),ones(p+r,1));
    se[1:p] = sqrt(diag(ioin));
  else;
    se = sqrt(diag(invpd(ion)));
  endif;

  format 8,3;
  fmt = {"*.*s" 8 8};
  t = {"*.*lf" 10 4};
  j = 1;
  do until j == 4;
    fmt = fmt|t;
    j = j + 1;
  endo;
  r = rows(s);
  names = name |"s0";
  t = se;
  t[p+1] = se[p+1]/sqrt(s0)/2;
  j = 1;
  do until j == r+ 1;
    t[p+j+1] = se[p+j+1]/sqrt(s[j])/2;
    names = names | "s"$+ftocv(j,1,0);
    j = j + 1;
  endo;

  format 8,4;

  /* transform from variance components to standard deviations */
  mask = zeros(rows(names),1)~ones(rows(names),rows(fmt)-1);
  ?;"Variable  Coefficient   SE      Est/SE ";
  ret = printfm(names~(a|sqrt(s0|s))~t~((a|sqrt(s0|s))./t),mask,fmt);?;
 "  Table reports standard deviations instead of variance component ests";
  ?;"Variance component estimates";
  "       Est        SE";
  ret = printfm(s~se[p+1:p+r],mask[p+1:p+r,2 3],fmt[2 3,.]);?;

  format 1,0;"  N iterations = ";;rows(iter);format 8,2;
  /* calc. likelihood */
  va = eigh(Pm);
  t = -(1/2)*( y'Pm*y - sumc(ln(selif(va,va.>0.001))) );
  format 8,3;"Log-likelihood = ";;t;

  output off;

RETP(ion,se);
ENDP;

PROC(3) = MEMN(a,s,r,z,q,x,y,nvec,respv,eff,s,test,lflag);
  output on;
  ?;"********* MEM algorithm for GLMM normal random effects *************";
  output off;
  if r /= rows(s);
    "s0 is wrong length";
  endif;

  /* setup indexes for random effects */
  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  /* initialize */
  n = rows(y);
  b = zeros(qs,1);
  p = cols(x);
  adj2 = zeros(r,1);
  nit = 0;
  iter = {};

  if lflag;
    locate 10,1;"iter    llk     a ... s  ";
  else;
    locate 10,1;"iter    a ... s  ";
  endif;
  /* start of MEM.  MEM cycles until converge = 1 */
  converge = 0;
  do until converge;

    olpv = pv;

    /* build matrix D and D inverse */
    I = eye(qs);
    d = {};
    j = 1;
    do until j == r + 1;
      t = seqa(qidx[j,1],1,q[j]);
      d = d | s[j]*ones(q[j],1);
      j = j + 1;
    endo;
    din = diagrv(I,1./d);

    /* compute conditional mean (m) and diagonal of conditional variance (w) */
    m = M_B(respv,(x*a + z*b));
    w = W_B(respv,m);

    /* Newton Raphson algorithm for computing b(p) */
    conve = 1;
    do until conve < 1E-5;

      vin = invpd(z'(w.*z) + din);
      adj1 = vin*(z'(y - m) - din*b);

      b = b + adj1;
      m = M_B(respv,x*a + z*b);
      w = W_B(respv,m);
      conve = sumc(abs(b - bi));
      bi = b;
    endo;
    if lflag;
      @ lk is Breslow and Clayton's (1993) approx. logL @
      if respv $== "P";
        lk = (x*a + z*b)'y - sumc(m) - sumc(lnfact(y)) - (1/2)*b'din*b
            - (1/2)*ln(det((d.*z')*(w.*z)+eye(qs)));
      else;
        lk = (x*a + z*b)'y + nvec'ln((nvec - m)./nvec) - (1/2)*b'din*b
           - (1/2)*ln(det((d.*z')*(w.*z)+eye(qs)));
      endif;
    endif;

    /* compute quantities for Laplace approximations */
    qvec = Q_B(respv,m);  @ second derivative of m wrt to eta @
    qzvin = qvec.*z*vin;

    /* compute Laplace approximation for X'm */
    adj2 = A_NXM(x,b,vin,qzvin,w,qvec);
    xm = x'm + adj2;

    /* M-step for alp */
    a = a + invpd(x'(x.*w))*(x'y - xm);

    /* compute Laplace approximation for bk'bk 's */
    adj2 = A_NT(b,vin,qzvin,qidx);

    /* M-step for s */
    j = 1;
    do until j == r + 1;
      t = seqa(qidx[j,1],1,q[j]);
      s[j] = (b[t]'b[t] + adj2[j])/q[j];
      j = j + 1;
    endo;

    /* pv is parameter vector */
    pv = s|a;
    converge = rows(pv) == sumc(abs(pv-olpv) .< test);
    nit = nit + 1;

    locate 11,1;
    format 4,0;nit;;format 7,3;
    if lflag;
      lk;;a';;s';
    else;
      a';;s';
    endif;
    st = st|t;  @ vector of approx. logL's@
    iter = iter|(a'~s');  @ save iterates @

  endo;
RETP(a,s,b);
ENDP;

PROC(3) = MEMLG(a,s,r,z,q,x,y,nvec,respv,eff,s,sn,test,lflag);
  output on;
  ?;"********* MEM algorithm for GLMM log-gamma effects *************";
  output off;

  if lflag;
    locate 10,1;"iter    llk     a ... s  ";
  else;
    locate 10,1;"iter    a ... s  ";
  endif;

  /* setup indexes for random effects */
  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  /* initialize */
  n = rows(y);
  p = cols(x);
  if r /= rows(s);
    "s0 is wrong length";
  endif;
  b = zeros(qs,1);
  adj2 = zeros(r,1);
  iter = {};

  converge = 0;
  do until converge;

    olpv = pv;

    m = M_B(respv,x*a + z*b);
    w = W_B(respv,m);

    /* build precursor (e) to Hessian matrix for computing b(p) */
    e = zeros(qs,1);
    i = 1;
    do until i == r + 1;
      t = seqa(qidx[i,1],1,q[i]);
      e[t] = s[i]*ones(q[i],1);
      i = i + 1;
    endo;

    /* Newton Raphson algorithm for computing b(p) */
    bi = b;
    I = eye(qs);
    conve = 1;
    do until conve < 1E-5;
      vin = invpd((z.*w)'z + diagrv(I,e.*exp(sn*b)));
      adj = vin*(z'(y - m) + sn*e.*(1 - exp(sn*b)));

      b = b + adj;
      m = M_B(respv,x*a + z*b);
      w = W_B(respv,m);
      conve = sumc(abs(b - bi));
      bi = b;
    endo;
    if lflag;
      if respv $== "P";
        lk = (x*a + z*b)'y - sumc(m) - sumc(lnfact(y))
            + q'(s.*ln(s) - ln(gamma(s))) + e'(sn*b - exp(sn*b))
            + (1/2)*ln(det(vin));
      else;
        lk = (x*a + z*b)'y + nvec'ln((nvec - m)./nvec)
            + q'(s.*ln(s) - ln(gamma(s))) + e'(sn*b - exp(sn*b))
            + (1/2)*ln(det(vin));
      endif;
    endif;

    /* compute quantities for Laplace approximations */
    qvec = Q_B(respv,m);
    qzvin = qvec.*z*vin;

    /* compute Laplace approximation for X'm */
    adj2 = A_LGXM(b,vin,qzvin,w,qvec,q,qidx);
    xm = x'm + adj2;

    /* M-step for alp */
    a = a + invpd(x'(x.*w))*(x'y - xm);

    /* compute Laplace approximation for s */
    adj2 = A_LGT(s,b,vin,qzvin,q,qidx);

    /* M-step for s */
    s = LG_E(s,b,adj2);

    /* parameter vector */
    pv = s|a;

    locate 11,1;
    format 4,0;nit;;format 7,3;
    if lflag;
      lk;;a';;1./s';
    else;
      a';;1./s';
    endif;

    converge = rows(pv) == sumc(abs(pv-olpv) .< test);
    nit = nit + 1;
    iter = iter|(a'~s');  @ save iterates @
  endo;
  output on;
  "Fixed effects estimates";format 8,4;a';
  "Variance component estimates";format 8,4;1./s';
  "Scale parameter (gamma) ";format 8,4;s';
  output off;
RETP(a,s,b);
ENDP;

PROC(2) = IOMEMN(a,b,s,respv,x,z,q,name,dname);
  @ computes approx. information matrix by differentiating the
     Laplace approximation of the marginal likelihood (see Steele [1996]) @
  output on;
  ?;"SE Computation for GLMM Model for normal random effects";
  $dname;

  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;

  p = cols(x);
  r = rows(qidx);

  m = M_B(respv,x*a + z*b);
  w = W_B(respv,m);
  qvec = Q_B(respv,m);
  qzvin = qvec.*z*vin;

  I = eye(qs);
  d = {};
  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    d = d | s[j]*ones(q[j],1);
    j = j + 1;
  endo;
  din = diagrv(I,1./d);
  vin = invpd(z'(w.*z) + din);

  ion = zeros(p+r,p+r);
  ion[1:p,1:p] = X'(W.*X)-X'(W.*Z*vin*Z'(W.*X));
  k = 1;
  do until k == r + 1;
    t = seqa(qidx[k,1],1,q[k]);
    Ik = zeros(qs,qs);
    Ik[t,t] = eye(q[k]);
    ion[1:p,p+k] = (s[k]^-2)*X'(W.*Z)*vin*Ik*b;
    ion[p+k,1:p] = ion[1:p,p+k]';

    ion[p+k,p+k] = (s[k]^-2)*q[k]/2 - (s[k]^-4)*b'Ik*vin*Ik*b;
    t = seqa(qidx[k,1],1,q[k]);
    Ik = zeros(qs,qs);
    Ik[t,t] = eye(q[k]);
    Dw = qzvin*Ik*b;
    if respv $== "P";
      U = m;
    else;
      U = m.*(1 - m./nvec).*(1 - 6*(m./nvec).*(1 - m./nvec));
    endif;
    e = (s[k]^-4)*tr(vin*(Ik - Z'(Dw.*Z))*vin*Ik)
       -(s[k]^-4)*tr(vin*(Ik - Z'(Dw.*Z))*vin*Z'(Dw.*Z))
       -(s[k]^-4)*tr(vin*Z'((((U.*z*vin*Ik*b).*Z)*vin*Ik*b).*Z))
       -(s[k]^-4)*tr(vin*(Z'((qzvin*(Ik - Z'(Dw.*Z))*vin*Ik*b).*Z)));
    ion[p+k,p+k] = ion[p+k,p+k] - e/2;

    i = k + 1;
    do until i == r + 1;
      Ij = zeros(qs,qs);
      t = seqa(qidx[i,1],1,q[i]);
      Ij[t,t] = eye(q[i]);
      ion[p+i,p+k] = (1/2)*(s[k]^-2)*(b'Ik*vin*Ij*b)*(s[i]^-2);

      ak = (s[k]^-2)*(s[i]^-2)*tr(vin*Ij*vin*(Ik - Z'(Dw.*Z)))
          - (s[k]^-2)*(s[i]^-2)*tr(vin*(Z'((qzvin*Ik*vin*Ij*b).*Z)));
      ion[p+i,p+k] = ion[p+i,p+k] - ak/2;
      ion[p+k,p+i] = ion[p+i,p+k];
      i = i + 1;
    endo;

    k = k + 1;
  endo;

  oldval = trapchk(1);
  trap 1,1;
  ioin = invpd(ion);
  trap oldval,1;
  if scalerr(ioin);
    "Info matrix singular";
    t = seqa(1,1,p);
    ioin = invpd(ion[t,t]);
    se = miss(ones(p+r,1),ones(p+r,1));
    se[1:p] = sqrt(diag(ioin));
  else;
    se = sqrt(diag(invpd(ion)));
  endif;

  format 8,3;
  fmt = {"*.*s" 8 8};
  t = {"*.*lf" 10 4};
  j = 1;
  do until j == 4;
    fmt = fmt|t;
    j = j + 1;
  endo;
  r = rows(s);
  names = name;
  t = se;
  j = 1;
  do until j == r+ 1;
    t[p+j] = se[p+j]/sqrt(s[j])/2;
    names = names | "s"$+ftocv(j,1,0);
    j = j + 1;
  endo;

  format 8,4;
  if respv $== "P";
    e = (x*a + z*b)'y - sumc(m) - (1/2)*b'din*b
        - (1/2)*ln(det((d.*z')*(w.*z)+eye(qs)));
  else;
    e = (x*a + z*b)'y + nvec'ln((nvec - m)./nvec) - (1/2)*b'din*b
        - (1/2)*ln(det((d.*z')*(w.*z)+eye(qs)));
  endif;
  ?;"Adjusted log-like = ";;e;

  /* transform from variance components to standard deviations */
  mask = zeros(rows(names),1)~ones(rows(names),rows(fmt)-1);
  ?;"Variable  Coefficient   SE      Est/SE ";
  ret = printfm(names~(a|sqrt(s))~t~((a|sqrt(s))./t),mask,fmt);?;
  "  Table reports standard deviations instead of variance component ests";
  format 1,0;"  N iterations = ";;rows(iter);format 10,4;
  ?;"Variance component estimates";
  "       Est        SE";
  ret = printfm(s~se[p+1:p+r],mask[p+1:p+r,2 3],fmt[2 3,.]);?;
  output off;
RETP(se,ion);
ENDP;

PROC(2) = IOMEMNF(a,b,s,respv,x,z,q,name,dname);
  @ computes approx. information matrix by differentiating the
     Laplace approximation of the marginal likelihood (see Steele [1996]) @
  output on;
  ?;"Short-cut SE Computation for GLMM Model";
  $dname;

  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;
  p = cols(x);
  r = rows(qidx);

  m = M_B(respv,x*a + z*b);
  w = W_B(respv,m);

  I = eye(qs);
  d = {};
  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    d = d | s[j]*ones(q[j],1);
    j = j + 1;
  endo;
  din = diagrv(I,1./d);
  vin = invpd(z'(w.*z) + din);

  ion = zeros(p+r,p+r);
  ion[1:p,1:p] = X'(W.*X)-X'(W.*Z*vin*Z'(W.*X));
  k = 1;
  do until k == r + 1;
    t = seqa(qidx[k,1],1,q[k]);
    Ik = zeros(qs,qs);
    Ik[t,t] = eye(q[k]);
    ion[1:p,p+k] = (s[k]^-2)*X'(W.*Z)*vin*Ik*b;
    ion[p+k,1:p] = ion[1:p,p+k]';

    ion[p+k,p+k] = (s[k]^-2)*q[k]/2 - (s[k]^-4)*b'Ik*vin*Ik*b;
    e = {};
    i = 1;
    do until i == r + 1;
      t = seqa(qidx[i,1],1,q[i]);
      e = e | b[t]*(i==j);
      i = i + 1;
    endo;
    t = seqa(qidx[k,1],1,q[k]);
    Ik = zeros(qs,qs);
    Ik[t,t] = eye(q[k]);

    i = k + 1;
    do until i == r + 1;
      Ij = zeros(qs,qs);
      t = seqa(qidx[i,1],1,q[i]);
      Ij[t,t] = eye(q[i]);
      ion[p+i,p+k] = (1/2)*(s[k]^-2)*(b'Ik*vin*Ij*b)*(s[i]^-2);
      ion[p+k,p+i] = ion[p+i,p+k];
      i = i + 1;
    endo;

    k = k + 1;
  endo;

  oldval = trapchk(1);
  trap 1,1;
  ioin = invpd(ion);
  trap oldval,1;
  if scalerr(ioin);
    "Info matrix singular";
    t = seqa(1,1,p);
    ioin = invpd(ion[t,t]);
    se = miss(ones(p+r,1),ones(p+r,1));
    se[1:p] = sqrt(diag(ioin));
  else;
    se = sqrt(diag(invpd(ion)));
  endif;

  format 8,3;
  fmt = {"*.*s" 8 8};
  t = {"*.*lf" 10 4};
  j = 1;
  do until j == 4;
    fmt = fmt|t;
    j = j + 1;
  endo;
  r = rows(s);
  names = name;
  t = se;
  j = 1;
  do until j == r+ 1;
    t[p+j] = se[p+j]/sqrt(s[j])/2;
    names = names | "s"$+ftocv(j,1,0);
    j = j + 1;
  endo;

  format 8,4;
  if respv $== "P";
    e = (x*a + z*b)'y - sumc(m) - (1/2)*b'din*b
        - (1/2)*ln(det((d.*z')*(w.*z)+eye(qs)));
  else;
    e = (x*a + z*b)'y + nvec'ln((nvec - m)./nvec) - (1/2)*b'din*b
        - (1/2)*ln(det((d.*z')*(w.*z)+eye(qs)));
  endif;
  ?;"Adjusted log-like = ";;e;

  /* transform from variance components to standard deviations */
  mask = zeros(rows(names),1)~ones(rows(names),rows(fmt)-1);
  ?;"Variable  Coefficient   SE      Est/SE ";
  ret = printfm(names~(a|sqrt(s))~t~((a|sqrt(s))./t),mask,fmt);?;
  "  Table reports standard deviations instead of variance component ests";
  format 1,0;"  N iterations = ";;rows(iter);format 10,4;
  ?;"Variance component estimates";
  "       Est        SE";
  ret = printfm(s~se[p+1:p+r],mask[p+1:p+r,2 3],fmt[2 3,.]);?;
  output off;
RETP(se);
ENDP;

PROC(2) = IOMEMLG(a,b,s,respv,x,z,q,name,dname);
  @ computes approx. information matrix by differentiating the
     Laplace approximation of the marginal likelihood (see Steele [1996]) @
  output on;
  ?;"SE Computation for GLMM Model with log-gamma random effects";
  $dname;

  qs = sumc(q);
  qidx = cumsumc(q);
  qidx = (trimr(1|qidx+1,0,1))~qidx;
  p = cols(x);
  r = rows(qidx);

  m = M_B(respv,x*a + z*b);
  w = W_B(respv,m);

  I = eye(qs);
  D = {};
  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    D = D | s[j]*ones(q[j],1);
    j = j + 1;
  endo;
  D = diagrv(I,D.*exp(sn*b));

  vin = invpd((z.*w)'z + D);

  ion = zeros(p+r,p+r);
  ion[1:p,1:p] = X'(W.*X) - X'(W.*Z*vin*Z'(W.*X));
  k = 1;
  do until k == r + 1;
    t = seqa(qidx[k,1],1,q[k]);
    Ik = zeros(qs,qs);
    Ik[t,t] = eye(q[k]);
    ion[1:p,p+k] = sn*X'(W.*Z)*vin*Ik*(ones(qs,1) - exp(sn*b));
    ion[p+k,1:p] = ion[1:p,p+k]';

    Dw = qzvin*Ik*(ones(qs,1) - exp(sn*b));
    U = qzvin*(exp(sn*b).*vin*Ik*(ones(qs,1) - exp(sn*b)));
    T = D.*vin.*(exp(sn*b).*vin*Ik*(ones(qs,1) - exp(sn*b)));

    Dg = -Ik*exp(sn*b);
    Dw = sn*qzvin*Ik*(1 - exp(sn*b));
    Dh = -D.*(vin*Ik.*(1 - exp(sn*b)));
    ion[p+k,p+k] = -q[k]*(1/s[k] - TRIGAM(s[k]))
          - (ones(qs,1) - exp(sn*b))'Ik*vin*Ik*(ones(qs,1) - exp(sn*b));

    i = k + 1;
    do until i == r + 1;
      Ij = zeros(qs,qs);
      t = seqa(qidx[i,1],1,q[i]);
      Ij[t,t] = eye(q[i]);
      ion[p+k,p+i] =
          - (ones(qs,1) - exp(sn*b))'Ik*vin*Ij*(ones(qs,1) - exp(sn*b));
      ion[p+i,p+k] = ion[p+k,p+i];

      ion[p+i,p+k] = ion[p+i,p+k] - ak/2;
      ion[p+k,p+i] = ion[p+i,p+k];
      i = i + 1;
    endo;

    k = k + 1;
  endo;

  oldval = trapchk(1);
  trap 1,1;
  ioin = invpd(ion);
  trap oldval,1;
  if scalerr(ioin);
    "Info matrix singular";
    t = seqa(1,1,p);
    ioin = invpd(ion[t,t]);
    se = miss(ones(p+r,1),ones(p+r,1));
    se[1:p] = sqrt(diag(ioin));
  else;
    se = sqrt(diag(invpd(ion)));
  endif;

  format 8,3;
  fmt = {"*.*s" 8 8};
  t = {"*.*lf" 10 4};
  j = 1;
  do until j == 4;
    fmt = fmt|t;
    j = j + 1;
  endo;
  r = rows(s);
  names = name;
  t = se;
  e = zeros(qs,1);
  j = 1;
  do until j == r+ 1;
    e[seqa(qidx[j,1],1,q[j])] = s[j]*ones(q[j],1);

    t[p+j] = se[p+j]*((s[j])^(-3/2))/2;
    names = names | "s"$+ftocv(j,1,0);
    j = j + 1;
  endo;

  format 8,4;
  if respv $== "P";
        e = (x*a + z*b)'y - sumc(m) - sumc(lnfact(y))
            + q'(s.*ln(s) - lnfact(s-1)) + e'(sn*b - exp(sn*b))
            + (1/2)*ln(det(vin));
  else;
        e = (x*a + z*b)'y + nvec'ln((nvec - m)./nvec)
            + q'(s.*ln(s) - lnfact(s-1)) + e'(sn*b - exp(sn*b))
            + (1/2)*ln(det(vin));
  endif;
  ?;"Adjusted log-like = ";;e;

  /* transform from variance components to standard deviations */
  mask = zeros(rows(names),1)~ones(rows(names),rows(fmt)-1);
  ?;"Variable  Coefficient   SE      Est/SE ";
  ret = printfm(names~(a|sqrt(1./s))~t~((a|sqrt(1./s))./t),mask,fmt);?;
  "  Table reports standard deviations instead of variance component ests";
  format 1,0;"  N iterations = ";;rows(iter);format 10,4;
  ?;"Variance component estimates";
  "       Est        SE";
  ret = printfm(s~se[p+1:p+r],mask[p+1:p+r,2 3],fmt[2 3,.]);?;
  output off;
RETP(se,ion);
ENDP;

PROC(1) = A_NXM(x,b,vin,qzvin,w,qvec);
  adj2 = zeros(cols(x),1);
  t = Z*vin*Z';
  Dw = qzvin*z'(w.*x);
  Dg = qvec.*X;
  j = 1;
  do until j == p + 1;
    adj2[j] = -(1/2)*tr(t.*(Dw[.,j]-Dg[.,j]));
    j = j + 1;
  endo;
RETP(adj2);
ENDP;

PROC(1) = A_LGXM(b,vin,qzvin,w,qvec,q,qidx);
  local e,t,i,j;
  Dh = {};
  r = rows(q);
  Dh = -exp(sn*b).*diag(vin);

  adj2 = zeros(p,1);
  t = z*vin*z';
  j = 1;
  do until j == p + 1;
    Dw = qzvin*z'(w.*x[.,j]);
    adj2[j] = -(1/2)*tr(t.*Dw) + (1/2)*tr(vin*z'(qvec.*x[.,j].*z))
                 + (1/2)*diag(vin)'(Dh.*(z'(w.*x[.,j])));
    j = j + 1;
  endo;
RETP(adj2);
ENDP;

PROC(1) = A_LGT(s,b,vin,qzvin,q,qidx);
  local i,j;
  adj2 = zeros(r,1);
  qs = sumc(q);
  r = rows(qidx);

  I = eye(qs);
  D = {};
  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    D = D | s[j]*ones(q[j],1);
    j = j + 1;
  endo;
  D = diagrv(I,D.*exp(sn*b));

  j = 1;
  do until j == r + 1;
    t = seqa(qidx[j,1],1,q[j]);
    Ij = zeros(qs,qs);
    Ij[t,t] = eye(q[j]);
    Dg = -Ij*exp(sn*b);
    Dw = sn*qzvin*Ij*(1 - exp(sn*b));
    Dh = -D.*(vin*Ij.*(1 - exp(sn*b)));
    adj2[j] = -(1/2)*tr(Z*vin*(Z.*Dw)') + (1/2)*tr(vin.*Dg) + (1/2)*tr(vin.*Dh);
    j = j + 1;
  endo;
RETP(adj2);
ENDP;

PROC(1) = M_B(respv,eta);
  if respv $== "P";
    m = exp(eta);
  else;
    m = nvec./(1 + exp(-eta));
  endif;
RETP(m);
ENDP;

PROC(1) = W_B(respv,m);
  if respv $== "P";
    w = m;
  else;
    w = m.*(1 - m./nvec);
  endif;
RETP(w);
ENDP;

PROC(1) = Q_B(respv,m);
  if respv $== "P";
    qvec = m;
  else;
    qvec = m.*(1 - m./nvec).*(1 - 2*m./nvec);
  endif;
RETP(qvec);
ENDP;

PROC(1) = MEMINIT(y,x,nvec,respv,sn);
  format /m1/rd 8,4;

  output on;
  if eff $== "LG";
    if sn == 1;
      "Log gamma random effects.";
    else;
      "Negative Log gamma random effects.";
    endif;
  else;
    "Normal random effects.";
  endif;
  if respv $== "P";
    "Conditional Poisson response variable. ";output off;
    a = PNEWT(y,x);
  else;
    "Conditional Binomial response variable. ";output off;
    a = BNEWT(y,x,nvec);
  endif;

RETP(a);
ENDP;

PROC(1) = A_NT(b,vin,qzvin,qidx);
  r = rows(qidx);
  adj2 = zeros(r,1);
  j = 1;
  do until j == r + 1;
    e = {};
    i = 1;
    do until i == r + 1;
      t = seqa(qidx[i,1],1,q[i]);
      e = e | 2*b[t]*(i==j);
      i = i + 1;
    endo;
    t = seqa(qidx[j,1],1,q[j]);
    Dw = qzvin*e;
    adj2[j] = -(1/2)*tr((Z.*Dw)'z*vin) + (1/2)*tr(2*vin[t,t]);
    j = j + 1;
  endo;
RETP(adj2);
ENDP;

PROC(1) = A_LGb(s,b,qzvin,veczvin);
  /* computes second order adjustment for approximating E[b|y] for
       log-gamma effects */

  local j,ii;
  adj1 = zeros(qs,1);
  Dh = -sn*exp(sn*b).*vin;
  ij = 0;
  j = 1;
  do until j == r + 1;
    ii = 1;
    do until ii == q[j] + 1;
      ij = ij + 1;
      adj1[ij] = -(1/2)*veczvin'vec(qzvin[.,ij].*z)
                + (1/2)*s[j]*diag(vin)'Dh[.,ij];
      ii = ii + 1;
    endo;
    j = j + 1;
  endo;
RETP(adj1);
ENDP;


PROC(1) = A_Nb(qzvin,veczvin);
  /* computes second order adjustment for approximating E[b|y] for
       normal random effects */
  adj1 = zeros(qs,1);
  ij = 0;
  j = 1;
  do until j == r + 1;
    ii = 1;
    do until ii == q[j] + 1;
      ij = ij + 1;
      adj1[ij] = -(1/2)*veczvin'vec(qzvin[.,ij].*z);
      ii = ii + 1;
    endo;
    j = j + 1;
  endo;
RETP(adj1);
ENDP;


PROC(1) = LG_E(s,b,adj);
  local e,i;
  i = 1;
  do until i == rows(q) + 1;
    e = seqa(qidx[i,1],1,q[i]);
    t = sumc(sn*b[e] - exp(sn*b[e])) + adj[i];
    conve = 1;
    do until conve < 1E-5;
      ol = s[i];
      s[i] = s[i]
         - (ln(s[i]) + 1 - DIGAM(s[i]) + t/q[i])/(1./s[i] - TRIGAM(s[i]));
      conve = abs(1/s[i] - 1/ol);
    endo;
    i = i + 1;
  endo;
RETP(s);
ENDP;

PROC(1) = BNEWT(y,x,nvec);
  converge = 0;
  a = invpd(x'x)*x'ln(((y+1)./(nvec-y+1)));
  m = nvec./(1+exp(-x*a));
  w = m.*(1 - m./nvec);
  do until converge;
    ol = a;
    a = a + invpd(x'(x.*w))*x'(y - m);
    m = nvec./(1+exp(-x*a));
    w = m.*(1 - m./nvec);
    converge = sumc(abs(a - ol)) < 0.0001;
  endo;
  RETP(a);
ENDP;

PROC(1) = PNEWT(y,x);
  converge = 0;
  a = invpd(x'x)*x'ln(y + (y .== 0));
  m = exp(x*a);
  do until converge;
    ol = a;
    a = a + invpd(x'(x.*m))*x'(y - m);
    m = exp(x*a);
    converge = abs(a - ol) < 0.0001;
  endo;
  RETP(a);
ENDP;

PROC(1) = TR(X);
  RETP(sumc(diag(x)));
ENDP;

