In jean997/sherlockAsh: CAUSE: Causal Analysis Using Summary Effect Estimates

knitr::opts_chunk$set(echo = TRUE)

Notes about code in src/ll_funcs_v7.cpp

ll_v7_ij

This function calculates $\log P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}, \sigma_{1}^2, \sigma_2^2, \rho]$

$$ P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}, \sigma_{1}^2, \sigma_2^2] = (1-q) \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma) \Sigma A(\gamma)^T + S_i(\rho) \right) + \ q \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma + \eta) \Sigma A(\gamma + \eta)^T + S_i(\rho) \right) \ $$

Where $$ A(x) \Sigma A(x)^T + S_i(\rho)= \begin{pmatrix} 1 & 0 \ x & 1 \end{pmatrix} \begin{pmatrix} \sigma_{1}^2 & 0 \ 0 & \sigma_2^2 \end{pmatrix} \begin{pmatrix} 1 & x \ 0 & 1 \end{pmatrix} + S_i(\rho) \ = \begin{pmatrix} 1 & 0 \ x & 1 \end{pmatrix} \begin{pmatrix} \sigma_1^2 & x\sigma_1^2 \ 0& \sigma_2^2 \end{pmatrix} + S_i(\rho) \ = \begin{pmatrix} \sigma_1^2 & x\sigma_1^2 \ x\sigma_1^2 & x^2 \sigma_1^2 + \sigma_2^2 \end{pmatrix} + S_i (\rho)\ = \begin{pmatrix} \sigma_1^2 + s_1^2 & x\sigma_1^2 + \rho s_1 s_2 \ x\sigma_1^2 + \rho s_1 s_2 & x^2 \sigma_1^2 + \sigma_2^2 + s_2^2 \end{pmatrix} $$

Specify the variance matrix for the first part of the equation

 double v11_1 = R_pow_di(sigma1,2) + R_pow_di(s1, 2) ;
 double v12_1 = g*R_pow_di(sigma1,2) +  rho*s1*s2 ;
 double v22_1 = R_pow_di(sigma2,2) + R_pow_di(g*sigma1,2) +  R_pow_di(s2, 2) ;

Here g is $\gamma$. To calulate $N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma) \Sigma A(\gamma)^T + S_i(\rho) \right)$ we use

$$ N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, \begin{pmatrix} v_{11} & v_{12} \ v_{21} & v_{22} \end{pmatrix} \right) = N(\hat{\beta}{1i}; 0, v{11} ) N(\hat{\beta}{2i}; \hat{\beta}{1i}v_{21}/v_{11}, v_{22}-v_{21}v_{12}/v_{11}) $$

Calculate conditional mean and variance for $\hat{\beta}{21}$ (note that $v{12} = v_{21}$ in this case).

  double m2_g1_1 = b1*v12_1/v11_1;
  double v2_g1_1 = v22_1 - R_pow_di(v12_1,2)/v11_1;

Calculate $\log \left( (1-q) \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma) \Sigma A(\gamma)^T + S_i(\rho) \right) \right)$

lik1_ij = log(1-q) + Rf_dnorm4(b1, 0, sqrt(v11_1), 1) +
    Rf_dnorm4(b2, m2_g1_1, sqrt(v2_g1_1), 1);

Calculate $\log \left( (q) \cdot N\left( \begin{pmatrix} \hat{\beta}{1i}\ \hat{\beta}{2i} \end{pmatrix}; \begin{pmatrix} 0 \ 0 \end{pmatrix}, A(\gamma + \eta) \Sigma A(\gamma + \eta)^T + S_i(\rho) \right) \right)$. Here gp is $\eta$ because $\eta$ was called $\gamma^\prime$ in an earlier version of notes.

  double v11_2 = R_pow_di(sigma1,2) + R_pow_di(s1, 2) ;
  double v12_2 = (g + gp)*R_pow_di(sigma1,2) +  rho*s1*s2 ;
  double v22_2 = R_pow_di(sigma2,2) + R_pow_di((g+gp)*sigma1,2) +  R_pow_di(s2, 2) ;
  double m2_g1_2 = b1*v12_2/v11_2;
  double v2_g1_2 = v22_2 - R_pow_di(v12_2,2)/v11_2;
  double lik2_ij = log(q) + Rf_dnorm4(b1, 0, sqrt(v11_2), 1) +
    Rf_dnorm4(b2, m2_g1_2, sqrt(v2_g1_2), 1);

Add them together

double lik_ij = Rf_logspace_add(lik1_ij, lik2_ij);
return lik_ij;

ll_v7_i

This function integrates over the distribution of $\sigma_1$ and $\sigma_2$ so we calculate

$$ P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}] = \sum_{k=1}^{K} \pi_k P[\hat{\beta}{1i}, \hat{\beta}{2i} \vert q, \gamma, \eta, \hat{s}{1i}, \hat{s}{2i}, \sigma_{1k}^2, \sigma_{2k}^2] $$

  int k = sigma1.length();
  std::vector<double> lik_i(k);
  for(size_t j = 0; j < k; j ++){
    lik_i[j] = log(pi[j]) +
      ll_v7_ij(rho, g, gp, q,
              sigma1[j], sigma2[j],
              b1,  b2, s1, s2);
  }
  double result = Rf_logspace_sum(&lik_i[0], k);
  return result;

ll_v7

Adds up over variants

  int p = beta_hat_2.size();
  NumericVector tll(p);
  RVector<double> ll(tll);

  parallel_for(blocked_range<size_t>(0, p),
               [&](const blocked_range<size_t>& r){
                 for(size_t i=r.begin(); i!=r.end(); ++i){

                   ll[i] = ll_v7_i(rho, g, gp, q,
                                   sigma1,sigma2,pi,
                                   beta_hat_1[i], beta_hat_2[i], seb1[i], seb2[i]);
                 }});
  double result = sum(tll);
  return result;

jean997/sherlockAsh documentation built on May 18, 2019, 11:45 p.m.