[DATALAD] Recorded changes

code/models/algo_siblings_adu.stan

@@ -102,7 +102,7 @@ parameters {
 
     // speech rates
     vector<lower=0>[n_classes] alpha_child_level; // variance across recordings for a given child
-    matrix<lower=0>[n_classes-1,n_corpora] alpha_corpus_level; // variance among children
+    array[2] matrix<lower=0>[n_classes-1,n_corpora] alpha_corpus_level; // variance among children
     matrix<lower=0>[n_classes-1,n_corpora] mu_corpus_level; // child-level average
     vector<lower=0>[n_classes-1] alpha_pop_level; // variance among corpora
     vector<lower=0>[n_classes] mu_pop_level; // population level averages
@@ -136,42 +136,47 @@ model {
     );
 
     vector [2] ll;
+    int distrib;
     for (c in 1:n_children) {
         // if there is sibling data
         if (child_siblings[c]>=0) {
+            distrib = child_siblings[c]>0?2:1;
+
             mu_child_level[c,1] ~ gamma(
-                alpha_corpus_level[1,corpus[c]],
-                (alpha_corpus_level[1,corpus[c]]/(mu_corpus_level[1,corpus[c]]*exp(
+                alpha_corpus_level[distrib,1,corpus[c]],
+                (alpha_corpus_level[distrib,1,corpus[c]]/(mu_corpus_level[1,corpus[c]]*exp(
                    child_siblings[c]>0?beta_sib_och:0 
                 )))
             );
             mu_child_level[c,2:] ~ gamma(
-                alpha_corpus_level[2:,corpus[c]],
-                (alpha_corpus_level[2:,corpus[c]]./mu_corpus_level[2:,corpus[c]]*exp(
+                alpha_corpus_level[distrib,2:,corpus[c]],
+                (alpha_corpus_level[distrib,2:,corpus[c]]./mu_corpus_level[2:,corpus[c]]*exp(
                    child_siblings[c]>0?beta_sib_adu:0 
                 ))
             );
         }
         // otherwise
         else {
+            // assuming no sibling
             ll[1] = log(p_sib)+gamma_lpdf(
-                mu_child_level[c,1] | alpha_corpus_level[1,corpus[c]], alpha_corpus_level[1,corpus[c]]/(mu_corpus_level[1,corpus[c]]*exp(beta_sib_och))
+                mu_child_level[c,1] | alpha_corpus_level[2,1,corpus[c]], alpha_corpus_level[2,1,corpus[c]]/(mu_corpus_level[1,corpus[c]]*exp(beta_sib_och))
             );
             ll[1] += gamma_lpdf(
-                mu_child_level[c,2] | alpha_corpus_level[2,corpus[c]], alpha_corpus_level[2,corpus[c]]/(mu_corpus_level[2,corpus[c]]*exp(beta_sib_adu))
+                mu_child_level[c,2] | alpha_corpus_level[2,2,corpus[c]], alpha_corpus_level[2,2,corpus[c]]/(mu_corpus_level[2,corpus[c]]*exp(beta_sib_adu))
             );
             ll[1] += gamma_lpdf(
-                mu_child_level[c,3] | alpha_corpus_level[3,corpus[c]], alpha_corpus_level[3,corpus[c]]/(mu_corpus_level[3,corpus[c]]*exp(beta_sib_adu))
+                mu_child_level[c,3] | alpha_corpus_level[2,3,corpus[c]], alpha_corpus_level[2,3,corpus[c]]/(mu_corpus_level[3,corpus[c]]*exp(beta_sib_adu))
             );
 
+            // assuming sibling
             ll[2] = log(1-p_sib)+gamma_lpdf(
-                mu_child_level[c,1] | alpha_corpus_level[1,corpus[c]], alpha_corpus_level[1,corpus[c]]/(mu_corpus_level[1,corpus[c]])
+                mu_child_level[c,1] | alpha_corpus_level[1,1,corpus[c]], alpha_corpus_level[1,1,corpus[c]]/(mu_corpus_level[1,corpus[c]])
             );
             ll[2] += gamma_lpdf(
-                mu_child_level[c,2] | alpha_corpus_level[2,corpus[c]], alpha_corpus_level[2,corpus[c]]/(mu_corpus_level[2,corpus[c]])
+                mu_child_level[c,2] | alpha_corpus_level[1,2,corpus[c]], alpha_corpus_level[1,2,corpus[c]]/(mu_corpus_level[2,corpus[c]])
             );
             ll[2] += gamma_lpdf(
-                mu_child_level[c,3] | alpha_corpus_level[3,corpus[c]], alpha_corpus_level[3,corpus[c]]/(mu_corpus_level[3,corpus[c]])
+                mu_child_level[c,3] | alpha_corpus_level[1,3,corpus[c]], alpha_corpus_level[1,3,corpus[c]]/(mu_corpus_level[3,corpus[c]])
             );
             target += log_sum_exp(ll);
         }
@@ -184,7 +189,8 @@ model {
     alpha_pop_level ~ gamma(8, 4); // sd = 0.35 x \mu
     alpha_pop ~ gamma(10, 10);
     for (i in 1:n_classes-1) {
-        alpha_corpus_level[i,:] ~ gamma(4, 4/alpha_pop[i]);
+        alpha_corpus_level[1,i,:] ~ gamma(4, 4/alpha_pop[i]);
+        alpha_corpus_level[2,i,:] ~ gamma(4, 4/alpha_pop[i]);
         mu_corpus_level[i,:] ~ gamma(alpha_pop_level[i],alpha_pop_level[i]/mu_pop_level[i+1]);
     }
 
@@ -205,16 +211,18 @@ model {
     }
 
     for (c in 1:n_speech_rate_children) {
+        distrib = child_siblings[c]>0?2:1;
+
         speech_rate_child_level[c,1] ~ gamma(
-            alpha_corpus_level[1,speech_rate_child_corpus[c]],
-            (alpha_corpus_level[1,speech_rate_child_corpus[c]]/(mu_corpus_level[1,speech_rate_child_corpus[c]]*exp(
+            alpha_corpus_level[distrib,1,speech_rate_child_corpus[c]],
+            (alpha_corpus_level[distrib,1,speech_rate_child_corpus[c]]/(mu_corpus_level[1,speech_rate_child_corpus[c]]*exp(
                 speech_rate_child_siblings[c]>0?beta_sib_och:0
             )))
         );
 
         speech_rate_child_level[c,2:] ~ gamma(
-            alpha_corpus_level[2:,speech_rate_child_corpus[c]],
-            (alpha_corpus_level[2:,speech_rate_child_corpus[c]]./(mu_corpus_level[2:,speech_rate_child_corpus[c]]*exp(
+            alpha_corpus_level[distrib,2:,speech_rate_child_corpus[c]],
+            (alpha_corpus_level[distrib,2:,speech_rate_child_corpus[c]]./(mu_corpus_level[2:,speech_rate_child_corpus[c]]*exp(
                 speech_rate_child_siblings[c]>0?beta_sib_adu:0
             )))
         );