disruption.stan 2.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899
  1. functions {
  2. vector z_scale(vector x) {
  3. return (x-mean(x))/sd(x);
  4. }
  5. }
  6. data {
  7. int<lower=1> N;
  8. int<lower=1> K;
  9. vector<lower=0>[N] soc_cap;
  10. vector<lower=0>[N] soc_div;
  11. vector<lower=0>[N] int_div;
  12. vector[N] res_soc_div;
  13. vector<lower=0>[N] age;
  14. vector<lower=0>[N] productivity;
  15. vector<lower=0>[N] productivity_solo;
  16. vector<lower=0>[N] m;
  17. matrix<lower=0,upper=1>[N,K] x;
  18. vector[N] stable;
  19. array [N] int<lower=0,upper=K-1> primary_research_area;
  20. }
  21. transformed data {
  22. vector[N] z_m = z_scale(m);
  23. vector[N] z_soc_cap = z_scale(soc_cap);
  24. vector[N] z_soc_div = z_scale(soc_div);
  25. vector[N] z_int_div = z_scale(int_div);
  26. vector[N] z_res_soc_div = z_scale(res_soc_div);
  27. vector[N] z_age = z_scale(age);
  28. vector[N] z_productivity = z_scale(productivity);
  29. vector[N] z_productivity_solo = z_scale(productivity_solo);
  30. }
  31. parameters {
  32. real beta_soc_cap;
  33. real beta_soc_div;
  34. real beta_int_div;
  35. real beta_stable;
  36. real beta_age;
  37. real beta_productivity;
  38. real beta_productivity_solo;
  39. vector[K] beta_x;
  40. real mu;
  41. real<lower=0> tau;
  42. real<lower=0> sigma;
  43. vector<lower=0,upper=1>[K] mu_x;
  44. vector<lower=1>[K] eta;
  45. real<lower=0,upper=1> mu_pop;
  46. real<lower=1> eta_pop;
  47. }
  48. model {
  49. vector[N] beta_research_area;
  50. for (k in 1:N) {
  51. beta_research_area[k] = beta_x[primary_research_area[k]+1]*tau;
  52. }
  53. beta_soc_cap ~ normal(0, 1);
  54. beta_soc_div ~ normal(0, 1);
  55. beta_int_div ~ normal(0, 1);
  56. beta_x ~ double_exponential(0, 1);
  57. beta_stable ~ normal(0, 1);
  58. beta_age ~ normal(0, 1);
  59. beta_productivity ~ normal(0, 1);
  60. beta_productivity_solo ~ normal(0, 1);
  61. mu ~ normal(0, 1);
  62. tau ~ exponential(1);
  63. sigma ~ exponential(1);
  64. z_m ~ normal(beta_soc_cap*z_soc_cap + beta_soc_div*z_res_soc_div + beta_int_div*z_int_div + beta_stable*stable + beta_age*z_age + beta_productivity*z_productivity + beta_productivity_solo*z_productivity_solo + beta_research_area + mu, sigma);
  65. eta ~ pareto(1, 1.5);
  66. mu_x ~ uniform(0, 1);
  67. eta_pop ~ pareto(1, 1.5);
  68. mu_pop ~ uniform(0, 1);
  69. for (k in 1:N) {
  70. //m[k] ~ beta_proportion(mu_x[primary_research_area[k]+1], eta[primary_research_area[k]+1]);
  71. }
  72. //m ~ beta_proportion(mu_pop, eta_pop);
  73. }
  74. generated quantities {
  75. real R2 = 0;
  76. {
  77. vector[N] beta_research_area;
  78. for (k in 1:N) {
  79. beta_research_area[k] = beta_x[primary_research_area[k]+1]*tau;
  80. }
  81. //vector[N] pred = inv_logit(beta_soc_cap*z_soc_cap + beta_soc_div*res_soc_div + beta_int_div*z_int_div + beta_stable*stable + beta_research_area + mu);
  82. vector[N] pred = beta_soc_cap*z_soc_cap + beta_soc_div*z_res_soc_div + beta_int_div*z_int_div + beta_stable*stable + beta_age*z_age + beta_productivity*z_productivity + beta_productivity_solo*z_productivity_solo + beta_research_area + mu;
  83. R2 = mean(square(z_m-pred))/variance(z_m);
  84. R2 = 1-R2;
  85. }
  86. }