target.cpp
Go to the documentation of this file.
1 
13 #include <limits>
14 #include <target/target.hpp>
15 #include <target/glm.hpp>
16 
17 namespace target {
18 
19  // typedef Target<T> B;
20 
33  template <typename T>
34  Target<T>::Target(const arma::Col<T> &y,
35  const arma::Mat<T> &a,
36  const arma::Mat<T> &x1,
37  const arma::Mat<T> &x2,
38  const arma::Mat<T> &x3,
39  const arma::Col<T> &parameter,
40  const arma::Col<T> &weights) : Target(y, a, x1, x2, x3, parameter) {
41  this->weights(weights);
42  }
43 
44 
50  template <typename T>
51  void Target<T>::update_par(const arma::Col<T> &parameter) {
52  unsigned pos = 0;
53  for (unsigned i=0; i < alpha.n_elem; i++) {
54  alpha[i] = parameter[i];
55  }
56  pos += alpha.n_elem;
57  for (unsigned i=0; i < beta.n_elem; i++) {
58  beta[i] = parameter[i+pos];
59  }
60  pos += beta.n_elem;
61  if (parameter.n_elem == pos+gamma.n_elem) {
62  for (unsigned i=0; i < gamma.n_elem; i++) {
63  gamma[i] = parameter[i+pos];
64  }
65  }
66  }
67 
68 
69  template <typename T>
70  Target<T>::Target(const arma::Col<T> &y,
71  const arma::Mat<T> &a,
72  const arma::Mat<T> &x1,
73  const arma::Mat<T> &x2,
74  const arma::Mat<T> &x3,
75  const arma::Col<T> &parameter) {
76  this->update_data(y, a, x1, x2, x3);
77  alpha = arma::Col<T>(x1.n_cols);
78  beta = arma::Col<T>(x2.n_cols);
79  gamma = arma::Col<T>(x3.n_cols);
80  update_par(parameter);
81  }
82 
83  template <typename T>
84  Target<T>::Target(const arma::Col<T> &y,
85  const arma::Mat<T> &a,
86  const arma::Mat<T> &x1,
87  const arma::Mat<T> &x2,
88  const arma::Col<T> &parameter,
89  const arma::Col<T> &weights) : Target(y, a, x1, x2, x2, parameter) {
90  this->weights(weights);
91  }
92 
93  template <typename T>
94  Target<T>::Target(const arma::Col<T> &y,
95  const arma::Mat<T> &a,
96  const arma::Mat<T> &x1,
97  const arma::Mat<T> &x2,
98  const arma::Col<T> &parameter) : Target(y, a, x1, x2, x2, parameter) {
99  }
100 
101 
102 
103  template <typename T>
104  void Target<T>::calculate(bool target, bool nuisance, bool propensity) {
105  if (target) // Target, linear predictor
106  this->target = Target<T>::X1()*Target<T>::alpha;
107  if (nuisance) // Nuisance, linear predictor
108  this->nuisance = Target<T>::X2()*Target<T>::beta;
109  if (propensity && Target<T>::gamma.n_elem > 0) { // Propensity
110  this->propensity = Target<T>::X3()*Target<T>::gamma;
111  this->propensity = target::expit(this->propensity);
112  }
113  }
114 
116 
124  template <typename T>
125  arma::Mat<T> TargetBinary<T>::pa() {
126  arma::Mat<T> res(this->Y().n_elem, 5);
127  res.col(1) = this->p(0);
128  res.col(2) = this->p(1);
129  res.col(0) = (1-this->A())%res.col(1) + this->A()%res.col(2);
130  res.col(3) = this->target;
131  res.col(4) = this->nuisance;
132  return(res);
133  }
134 
141  template <typename T>
142  arma::Col<T> TargetBinary<T>::loglik(bool indiv) {
143  arma::Col<T> phat = TargetBinary<T>::p(0) %
145  arma::Col<T> logl = Target<T>::Y() % log(phat) +
146  (1-Target<T>::Y()) % log(1-phat);
147  logl %= Target<T>::weights();
148  if (indiv) return(logl);
149  return(sum(logl, 0));
150  }
151 
152 
169  template <typename T>
170  arma::Mat<T> TargetBinary<T>::est(arma::Col<T> alpha,
171  const arma::Col<T> &propensity) {
172  arma::Col<T> p0 = this->p(0);
173  for (unsigned i=0; i < alpha.n_elem; i++) this->alpha[i] = alpha[i];
174  calculate(true, false, false);
175  // arma::Col<T> a = this->A().col(0);
176  arma::Col<T> H = this->H();
177  arma::Col<T> S = (this->A()-propensity)%(H-p0);
178  S %= this->weights();
179  arma::Mat<T> U(S.n_elem, alpha.n_elem);
180  for (unsigned i=0; i < alpha.n_elem; i++) {
181  U.col(i) = S % this->X1().col(i);
182  }
183  return(U);
184  }
185 
186  template <typename T>
187  arma::Mat<T> TargetBinary<T>::est(arma::Col<T> alpha) {
188  calculate(false, false, true); // Calculate propensity weights
189  return( sum(TargetBinary<T>::est(alpha, this->propensity), 0) );
190  }
191 
192 
199  template <typename T>
200  arma::Mat<T> TargetBinary<T>::score(bool indiv) {
201  arma::Col<T> phat = this->p(0) % (1-this->A().col(0)) +
202  this->p(1) % this->A().col(0);
203  arma::Mat<T> dp_dlp = dp();
204  arma::Col<T> S = (this->Y()-phat) / (phat%(1-phat));
205  S %= this->weights();
206  for (unsigned i=0; i < dp_dlp.n_cols; i++)
207  dp_dlp.col(i) %= S;
208  if (indiv) return(dp_dlp);
209  return(sum(dp_dlp, 0));
210  }
211 
212 
213  template <typename T>
214  void TargetBinary<T>::calculate(bool target,
215  bool nuisance,
216  bool propensity) {
217  Target<T>::calculate(target, nuisance, propensity);
218  if (nuisance)
219  this->nuisance = exp(this->nuisance); // Odds-product
220  }
221 
223 
224  template <typename T>
225  RR<T>::RR(const arma::Col<T> &y,
226  const arma::Mat<T> &a,
227  const arma::Mat<T> &x1,
228  const arma::Mat<T> &x2,
229  const arma::Mat<T> &x3,
230  const arma::Col<T> &parameter,
231  const arma::Col<T> &weights) :
232  TargetBinary<T>(y, a, x1, x2, x3, parameter, weights) {
233  RR::calculate();
234  }
235 
236  template <typename T>
237  RD<T>::RD(const arma::Col<T> &y,
238  const arma::Mat<T> &a,
239  const arma::Mat<T> &x1,
240  const arma::Mat<T> &x2,
241  const arma::Mat<T> &x3,
242  const arma::Col<T> &parameter,
243  const arma::Col<T> &weights) :
244  TargetBinary<T>(y, a, x1, x2, x3, parameter, weights) {
245  RD::calculate();
246  }
247 
248  template <typename T>
249  void RD<T>::calculate(bool target, bool nuisance, bool propensity) {
250  TargetBinary<T>::calculate(target, nuisance, propensity);
251  if (target)
252  this->target = tanh(this->target); // Relative risk
253  if (target || nuisance)
254  this->pr = rd2prob(this->target, this->nuisance);
255  }
256 
257  template <typename T>
258  void RR<T>::calculate(bool target, bool nuisance, bool propensity) {
259  TargetBinary<T>::calculate(target, nuisance, propensity);
260  if (target)
261  this->target = exp(this->target); // Relative risk
262  if (target || nuisance)
263  this->pr = rr2prob(this->target, this->nuisance);
264  }
265 
266  template <typename T>
267  arma::Mat<T> RD<T>::dp() {
268  arma::Col<T> a = RD<T>::A().col(0);
269  arma::Col<T> s0 = RD<T>::p(0) % (1 - RD<T>::p(0));
270  arma::Col<T> s1 = RD<T>::p(1) % (1 - RD<T>::p(1));
271  arma::Col<T> dp0_rd = -s0/(s0 + s1);
272  arma::Col<T> d0 = (dp0_rd + a) % (1 - rd()%rd());
273  arma::Mat<T> d_alpha = RD<T>::X1();
274  for (unsigned i=0; i < d_alpha.n_cols; i++) d_alpha.col(i) %= d0;
275  arma::Col<T> dp0_op = s0 % s1/(s0 + s1);
276  arma::Mat<T> d_beta = RD<T>::X2();
277  for (unsigned i=0; i < d_beta.n_cols; i++) d_beta.col(i) %= dp0_op;
278  return ( join_rows(d_alpha, d_beta) );
279  }
280 
281  template <typename T>
282  arma::Mat<T> RR<T>::dp() {
283  arma::Col<T> a = RR<T>::A().col(0);
284  arma::Col<T> s = (1 - RR<T>::p(0)) + (1 - RR<T>::p(1));
285  arma::Col<T> tmp = 1 - a + a%rr();
286  arma::Col<T> dp_rr = -RR<T>::p(0)/RR<T>::p(1) %
287  RR<T>::p(0)%(1-RR<T>::p(0)) / s % tmp + a % RR<T>::p(0);
288  dp_rr %= rr();
289  arma::Mat<T> d_alpha = RR<T>::X1();
290  for (unsigned i=0; i < d_alpha.n_cols; i++) d_alpha.col(i) %= dp_rr;
291  s = (1-RR<T>::p(1))%RR<T>::p(1) + (1-RR<T>::p(0))%RR<T>::p(1);
292  arma::Col<T> dp_op = pow((1-RR<T>::p(0))%(1-RR<T>::p(1)), 2) / s % tmp;
293  arma::Mat<T> d_beta = RR<T>::X2();
294  dp_op %= op();
295  for (unsigned i=0; i < d_beta.n_cols; i++) d_beta.col(i) %= dp_op;
296  return ( join_rows(d_alpha, d_beta) );
297  }
298 
306  template <typename T>
307  arma::Mat<T> rd2prob(const arma::Col<T> &rd, const arma::Col<T> &op) {
308  arma::Col<T> a = op-1;
309  arma::Col<T> b = -op%(2-rd)-rd;
310  arma::Col<T> p0 = (-b - sqrt(b%b - 4*op%(1.0 - rd)%a)) / (2*a);
311  for (unsigned i=0; i < p0.size(); i++) {
312  if (std::abs(op[i]-1.0) < 1e-16) p0[i] = 0.5*(1.0 - rd[i]);
313  }
314  arma::Mat<T> res(rd.n_elem, 2);
315  res.col(0) = p0;
316  res.col(1) = p0+rd;
317  return res;
318  }
319 
327  template <typename T>
328  arma::Mat<T> rr2prob(const arma::Col<T> &rr, const arma::Col<T> &op) {
329  arma::Mat<T> res(rr.n_elem, 2);
330  arma::Col<T> a = rr%(1.0 - op);
331  arma::Col<T> b = op%(1.0 + rr);
332  arma::Col<T> p0 = (-b + sqrt(b%b + 4*a%op))/(2*a);
333  for (unsigned i=0; i < p0.size(); i++) {
334  if (std::abs(op[i]-1.0) < 1e-16) p0[i] = (T)1/(1.0 + rr[i]);
335  }
336  res.col(0) = p0;
337  res.col(1) = p0%rr;
338  return res;
339  }
340 
342 
343 
357  ACE::ACE(const arma::cx_vec &y,
358  const arma::cx_mat &a,
359  const arma::cx_mat &x2,
360  const arma::cx_mat &x3,
361  const arma::cx_vec &parameter,
362  const arma::cx_vec &weights,
363  bool bin) :
364  Target<cx_dbl>(y, a, arma::cx_mat(1, 1), x2, x3, parameter, weights) {
365  this->binary = bin;
366  ACE::calculate(true, true);
367  }
368 
369  ACE::ACE(const arma::vec &y,
370  const arma::vec &a,
371  const arma::mat &x2,
372  const arma::mat &x3,
373  const arma::vec &parameter,
374  const arma::vec &weights,
375  bool bin) :
376  ACE(arma::conv_to<arma::cx_vec>::from(y),
377  arma::conv_to<arma::cx_vec>::from(a),
378  arma::conv_to<arma::cx_mat>::from(x2),
379  arma::conv_to<arma::cx_mat>::from(x3),
380  arma::conv_to<arma::cx_vec>::from(parameter),
381  arma::conv_to<arma::cx_vec>::from(weights),
382  bin) { }
383 
384  void ACE::calculate(bool target, bool nuisance, bool propensity) {
385  Target<cx_dbl>::calculate(false, nuisance, propensity);
386  if (nuisance && this->binary) {
387  this->nuisance = target::expit(this->nuisance); // outcome regression
388  }
389  }
390 
391  arma::cx_mat ACE::est(bool indiv, const cx_dbl &value) {
392  arma::cx_vec a1 = A().as_col();
393  for (unsigned i=0; i < a1.n_elem; i++) a1[i] = (a1[i] == value);
394  arma::cx_vec U = Y()%a1/(propensity) -
395  (a1-propensity)%nuisance/(propensity) - alpha[0];
396  U %= weights();
397  if (indiv) return( std::move(U) );
398  return(sum(U, 0));
399  }
400 
401  arma::cx_mat ACE::est(arma::cx_vec par, bool indiv, const cx_dbl &value) {
403  ACE::calculate();
404  return(ACE::est(indiv, value));
405  }
406 
407  void ACE::update_par(arma::cx_vec par) {
408  this->Target<cx_dbl>::update_par(par);
409  }
410 
411  void ACE::update_par(arma::vec par) {
412  arma::cx_vec parc = arma::conv_to<arma::cx_vec>::from(par);
413  this->Target<cx_dbl>::update_par(parc);
414  }
415 
423  arma::mat ACE::deriv(const cx_dbl &value) {
424  unsigned n1 = this->alpha.n_elem;
425  unsigned n2 = this->beta.n_elem;
426  unsigned n3 = this->gamma.n_elem;
427  unsigned p = n1+n2+n3;
428  arma::vec res(p);
429  // int n = Y().n_elem;
430  res[0] = -sum(real(weights()));
431  double eps = std::numeric_limits<float>::denorm_min();
432  cx_dbl delta(0, eps);
433  unsigned pos = 1;
434  cx_dbl tmp;
435 
436  for (unsigned i=0; i < n2; i++) {
437  tmp = this->beta[i];
438  this->beta[i] += delta;
439  ACE::calculate(true, true, false);
440  this->beta[i] = tmp;
441  res[i+pos] = std::imag(ACE::est(false, value)[0])/eps;
442  }
443  pos += n2;
444  for (unsigned i=0; i < n3; i++) {
445  tmp = this->gamma[i];
446  this->gamma[i] += delta;
447  ACE::calculate(i == 0, true);
448  this->gamma[i] = tmp;
449  res[i+pos] = std::imag(ACE::est(false, value)[0])/eps;
450  }
451  return( std::move(res) );
452  }
453 
454 
456 
457 
458  // specific template instantiations:
459  template class RR<double>;
460  template class RR<cx_dbl>;
461  template class RD<double>;
462  template class RD<cx_dbl>;
463  template class Target<double>;
464  template class Target<cx_dbl>;
465  template class TargetBinary<double>;
466  template class TargetBinary<cx_dbl>;
467 
468 
469 
470 } // namespace target
arma::mat deriv(const cx_dbl &value=1)
deriv Calculate derivative of estimating function for AIPW estimator
Definition: target.cpp:423
Classes for targeted inference models.
void update_par(const arma::Col< T > &parameter)
update_par -
Definition: target.cpp:51
arma::Mat< T > rr2prob(const arma::Col< T > &rd, const arma::Col< T > &op)
rr2prob - Computes risk probabilities given exposure 0 or 1.
Definition: target.cpp:328
arma::Mat< T > rd2prob(const arma::Col< T > &rd, const arma::Col< T > &op)
rd2prob - Computes risk probabilities given exposure 0 or 1.
Definition: target.cpp:307
Utility functions for Generalized Linear Models.
ACE(const arma::cx_vec &y, const arma::cx_mat &a, const arma::cx_mat &x2, const arma::cx_mat &x3, const arma::cx_vec &parameter, const arma::cx_vec &weights, bool binary=true)
ACE AIPW for Average Causal Effects.
Definition: target.cpp:357
Target()
Constructor.
Definition: target.hpp:38