Multivariate Normal Model for Heritability of Drift

 

Model for multivariate drift in allele frequencies away from 0.5.  This model takes account of the correlation in drift away from 0.5 by modelling the drift as a multivariate normal distribution on the logit of the allele frequencies.

 

We perform a hierarchical model because of the error structure.  The error is in two parts - teh first due to errors with the PCR reaction, the second due to errors with the gel runs.  These errors are hierarchical, so any classical analysis would have to use a split-plot type design - this would be very difficult as we have a lot of missing data - for most twin pairs we have either a second gel run or a second PCR and gel run, in very few do we have both.  A Bayesian analysis is easier with this much missing data.

 

Model

{

    # lp[twinpair][celltype][whichtwin][replicate]

    # replicate 2 is a repeated run

    # replicate 3 is repeated PCR and run

    # whichtwin = 1, 2 for first or second twin

    for (i in 1:N) {

             for (j in 1:2) {

                for (k in 1:2) {

                   for (m in 1:2) {

                   # both gel runs (if both occurred)

                   # the observed lp is normal mean lpp var 1/taue[1]

                      lp[i,j,k,m] ~ dnorm(lpp[i,j,k,1],taue[1]);

                   }

                   #  if we have a second PCR run then var 1/taue[1]

                   # lpp are the (missing) frequencies following PCR

                   lp[i,j,k,3] ~dnorm(lpp[i,j,k,2],taue[1]);

                   for (m in 1:2) {

                       lpp[i,j,k,m] ~ dnorm(lpm[i,j,k],taue[2]);

                   } 

              }

              lpm[i,j,1:2]~dmnorm(mu[j,],tau[j,code[i],,])

              }

      }

      # we have tau[celltype][twintype][twin1][twin2]

      # twintype==1 is for monozygotes;  twintype==2 dizygotes

      # where celltype==1  neutrophils;  celltype==2 tcells

      for (i in 1:2) {

         for (j in 1:2) {

             tau[i,j,1:2,1:2]~dwish(R[,],2);

             for (k in 1:2) {

                 for (m in 1:2) {

                     Sigma[i,j,k,m]<-inverse(tau[i,j,,],k,m);

                 }

                   #  F[i,j,k,k] <- Sigma[i,j,k,k]/(16*4);

             }

             #cov covariance[celltype][twintype]

             cor[i,j]<-Sigma[i,j,1,2]/(sqrt(Sigma[i,j,1,1]*Sigma[i,j,2,2]))

         }

         taue[i]~dgamma(0.001, 0.001)

         sigmae[i] <- 1/taue[i];

         h2[i] <- 2.*(cor[i,1]-cor[i,2]);

      }

}

 

Initial Values

list(taue=c(100,100),

lpm=structure(.Data=c(

-0.867500567704723, 0.0582689081239758, -0.245900538436826,

0.173953307123438, 1.10856261952128, -0.136965855073157, -0.908818717035454,

-0.266573109241546, -0.0565703514883942, -0.527632742082372,

0.548121408509687, 0.246860077931526, -2.39689577246529, -1.42295834549148,

0.0769610411361286, 0.371563556432483, -0.39452516806983, -0.359536176219765,

-1.27296567581289, -0.336872316642553, 1.02604159583327, 0.512823626428664,

0.536493370514568, 0.207014169384326, 0.444685821261446, 0.553885113226438,

0.0478373294141601, -0.0356271776431512, -0.221894331913778,

-0.104250021373799, 0.357674444271816, 0.500775287912489, -0.84863208340034,

-0.12897038129696, 0.916290731874155, 0.756121979721333, 0.559615787935423,

0.542324290825362, 0.165514438477574, 0.593326845277734, 0.405465108108164,

0.625938430866495, -0.294371060602578, -0.811930716549912, 0.165514438477574,

-0.207024169434327, 0.678033542749898, 0.62057648772511, 0.371563556432483,

0.65232518603969, -2.48891467118554, -0.638658995275876, 1.48160454092422,

-0.59783700075562, 0.29266961396282, 0.444685821261446, 0.285178942233662,

-0.161343150408763, 0.190620359608650, -0.0910193983871684, -1.57021719928082,

-1.86433016206289, -0.894040122939335, -0.449416995637347, 0.457424847038875,

-0.0387408283164307, -0.978166135592242, -0.918793862092274,

0.131028262406404, 0.378436435720245, -0.625488532086131, -0.49593701127224,

0.139761942375158, 0.148420005118273, 0.576613364303993, 0.806475865866949,

0.751416088683921, 0.165514438477574, 1.26129787094521, 0.444685821261446,

1.62136648329937, 1.44691898293633, 3.89995042419387, 4.0517849478033,

-0.98886142470899, -0.574475650842447, -0.644357016390513, -0.263965545834465,

0.357674444271816, 0.463734016232141, -0.0534007767271152, 0.0861776962410526

,-0.352398387171472, 0.095310179804325, -0.341082849178896,

0.173953307123438, 0.451075619360216, -0.424647927524938, -0.770028224895903,

-0.596020469829223, 1.85941811770187, 1.45861502269952, 0.712949807856125,

0.285178942233662, -2.22562405185792, -1.48722027970985, 1.41098697371026,

1.41098697371026, 1.83098018238134, 1.79840401194672, -1.49610922712710,

-0.262664309476493, 0.371563556432483, 0.542324290825362, 0.476234178996372,

0.207014169384326, 0.828551817566148, 0.779324876800998, -0.215671536475509,

-0.147340587898709, -0.3147107448397, -0.293029678778376, 0.887891257352457,

0.78845736036427, -0.116533816255952, 0.122217632724249, -3.10109278921182,

-2.76462055259060, 0.553885113226438, 0.506817602368451, -0.0212236364516268,

0.262364264467491, -0.223143551314210, -0.0408219945202552, -0.0910193983871684,

0.0295588022415442, -0.0704224642965458, -0.391562202939173,

0.565313809050061, 0.329303747142600, 0.451075619360216, 0.29266961396282,

-3.93222571274567, -3.29683736633791, 0.0198026272961797, 0.978326122793608,

2.13060982825424, 0.00995033085316831, 0.444685821261446, -0.105360515657826,

0.173953307123438, -0.00200200267067319, -1.73727128394399, -1.75446368448436,

-0.339677367570161, -0.400477566597125, 0.0676586484738147, 0.0861776962410526,

0.246860077931526, -0.173163619009189, 0.231111720963387, 0.048790164169432,

0.122217632724249, 0.122217632724249, -0.504181081047322, -0.211956361923646,

0.76080582903376, 0.633927820899974, -0.136965855073157, 0.0295588022415442,

1.41585316336144, 0.182321556793954, -0.211956361923646, -0.136965855073157,

0.683096844706444, 0.62057648772511, 1.97408102602201, 1.60542989103656,

0.683096844706444, 0.405465108108164, -0.00904074465214918, -0.415515443961666,

-0.805196684368568, -0.683196849706777,0.0,0.0,0.0,0.0), .Dim = c(47,2,2))

 

)

 

Data

list(

 N=47,

mu=structure(.Data=c(0,0,0,0), .Dim=c(2,2) ),

 

code = c(

1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

1, 1, 1, 1, 1, 1, 1, 1, 1, 1,

1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2,

2, 2, 2, 2, 2, 2, 2, 2, 2, 2,

2, 2, 2, 2, 2, 2),

 

R=structure(.Data=c(1,0,0,1),.Dim=c(2,2)),

 

lp = structure(.Data = c(-0.867500567704723, -0.772190387900398,

NA, -0.352398387171472, -0.398986142010455, NA, 0.0582689081239758,

0, NA, 0.095310179804325, 0.0487901641694320, NA, -0.245900538436826,

NA, -0.145025772050258, -0.341082849178896, NA, -0.302457358033935,

0.173953307123438, NA, 0.215111379616945, 0.173953307123438,

NA, 0.2390169004705, 1.10856261952128, 0.982078472412158, NA,

0.451075619360217, 0.350656871613169, NA, -0.136965855073157,

-0.126697653045958, NA, -0.424647927524938, -0.409473129505703,

NA, -0.908818717035454, -1.02722229258144, NA, -0.770028224895903,

-0.911303190363116, NA, -0.266573109241546, -0.55512588266257,

NA, -0.596020469829223, -0.703197516413447, NA, -0.0565703514883944,

0.104360015324243, NA, 1.85941811770187, 1.90805992492422, NA,

-0.527632742082372, -0.49429632181478, NA, 1.45861502269952,

1.49514876603197, NA, 0.548121408509688, 0.548121408509688, NA,

0.712949807856125, 0.708035793053696, NA, 0.246860077931526,

0.207014169384326, NA, 0.285178942233662, 0.285178942233662,

NA, -2.39689577246529, -2.25379492882461, NA, -2.22562405185792,

-2.13707065451647, NA, -1.42295834549148, -1.59454929994035,

NA, -1.48722027970985, -1.70925824771631, NA, 0.0769610411361284,

NA, 0.0582689081239758, 1.41098697371026, NA, 1.28093384546206,

0.371563556432483, NA, 0.277631736598280, 1.41098697371026, NA,

1.34807314829969, -0.39452516806983, -0.373966441048793, -0.369615455214467,

1.83098018238134, 1.54968790802833, 1.54968790802833, -0.359536176219765,

-0.368169323364468, -0.373966441048793, 1.79840401194672, 1.61143591509677,

1.61740608208328, -1.27296567581289, NA, -1.27296567581289, -1.49610922712710,

NA, -1.56542102701733, -0.336872316642553, NA, -0.436955775199535,

-0.262664309476493, NA, -0.355247391947547, 1.02604159583327,

NA, 1.02961941718116, 0.371563556432483, NA, 0.29266961396282,

0.512823626428664, NA, 0.476234178996372, 0.542324290825362,

NA, 0.542324290825362, 0.536493370514568, 0.553885113226438,

NA, 0.476234178996372, 0.457424847038875, NA, 0.207014169384326,

0.139761942375159, NA, 0.207014169384326, 0.148420005118273,

NA, 0.444685821261446, NA, NA, 0.828551817566148, NA, NA, 0.553885113226438,

NA, NA, 0.779324876800998, NA, NA, 0.0478373294141601, NA, NA,

-0.215671536475509, NA, NA, -0.0356271776431512, NA, NA, -0.147340587898709,

NA, NA, -0.221894331913778, NA, NA, -0.3147107448397, NA, NA,

-0.104250021373799, NA, NA, -0.293029678778376, NA, NA, 0.357674444271816,

0.512823626428664, NA, 0.887891257352457, 0.904218150639886,

NA, 0.500775287912489, 0.476234178996372, NA, 0.78845736036427,

0.756121979721334, NA, -0.84863208340034, NA, -0.911303190363116,

-0.116533816255952, NA, -0.123298216344494, -0.12897038129696,

NA, -0.219400565035375, 0.122217632724249, NA, 0.123985979780991,

0.916290731874155, 0.871293365943419, 0.727548607277278, -3.10109278921182,

-3.10109278921182, -3.2188758248682, 0.756121979721334, 0.966983846189673,

0.792992515529661, -2.76462055259060, -2.70306265959117, -2.86998106824843,

0.559615787935423, NA, NA, 0.553885113226438, NA, NA, 0.542324290825362,

NA, NA, 0.506817602368452, NA, NA, 0.165514438477573, NA, NA,

-0.0212236364516267, NA, NA, 0.593326845277734, NA, NA, 0.262364264467491,

NA, NA, 0.405465108108164, NA, 0.314810739840034, -0.223143551314210,

NA, -0.466808738349216, 0.625938430866495, NA, 0.392042087776024,

-0.0408219945202552, NA, -0.226900600191922, -0.294371060602578,

-0.365283318475333, -0.400477566597125, -0.0910193983871685,

-0.156653810045377, -0.182721636815294, -0.811930716549912, -0.811930716549912,

-0.814185508937001, 0.0295588022415444, 0.0582689081239758, -0.187535123846842,

0.165514438477573, NA, 0.0861776962410524, -0.0704224642965458,

NA, -0.0779615414697118, -0.207024169434327, NA, -0.136965855073157,

-0.391562202939173, NA, -0.207024169434327, 0.678033542749897,

NA, NA, 0.56531380905006, NA, NA, 0.62057648772511, NA, NA, 0.329303747142600,

NA, NA, 0.371563556432483, NA, 0.371563556432483, 0.451075619360217,

NA, 0.457424847038875, 0.65232518603969, NA, 0.494696241836107,

0.29266961396282, NA, 0.506817602368452, -2.48891467118554, -2.44184716032755,

NA, -3.93222571274567, -4.25451331437492, NA, -0.638658995275876,

-0.653926467406664, NA, -3.29683736633791, -3.36390159691846,

NA, 1.48160454092422, NA, 1.40044497500497, 0.0198026272961797,

NA, -0.0888312137066157, -0.59783700075562, NA, -0.59783700075562,

0.978326122793608, NA, 0.667829372575655, 0.29266961396282, NA,

0.246860077931526, 2.13060982825424, NA, 1.70837786028900, 0.444685821261446,

NA, 0.494696241836107, 0.0099503308531681, NA, -0.407968238326283,

-0.327116141697188, NA, NA, -0.269187489815617, NA, NA, -0.311974765020826,

NA, NA, -0.338273858567841, NA, NA, 0.285178942233662, NA, NA,

0.444685821261446, NA, NA, -0.161343150408763, NA, NA, -0.105360515657826,

NA, NA, 0.190620359608650, NA, NA, 0.173953307123438, NA, NA,

-0.0910193983871685, NA, NA, -0.00200200267067308, NA, NA, -1.57021719928082,

NA, NA, -1.73727128394399, NA, NA, -1.86433016206289, NA, NA,

-1.75446368448436, NA, NA, -0.894040122939335, NA, -0.881889305156823,

-0.339677367570161, NA, -0.515838165589535, -0.449416995637347,

NA, -0.387134151423441, -0.400477566597125, NA, -0.375420986759788,

0.457424847038875, NA, 0.500775287912489, 0.0676586484738149,

NA, 0.104360015324243, -0.0387408283164306, NA, -0.0110609473594249,

0.0861776962410524, NA, 0.165514438477573, -0.978166135592243,

-0.867500567704723, NA, 0.246860077931526, 0.215111379616945,

NA, -0.918793862092274, -0.711311151187616, NA, -0.173163619009189,

-0.20089294237939, NA, 0.131028262406404, 0.476234178996372,

NA, 0.231111720963387, 0.190620359608650, NA, 0.378436435720245,

0.536493370514568, NA, 0.0487901641694320, 0.122217632724249,

NA, -0.62548853208613, NA, NA, 0.122217632724249, NA, NA, -0.49593701127224,

NA, NA, 0.122217632724249, NA, NA, 0.139761942375159, NA, NA,

-0.504181081047322, NA, NA, 0.148420005118273, NA, NA, -0.211956361923645,

NA, NA, 0.576613364303994, NA, 0.438254930931155, 0.76080582903376,

NA, 1.07158361628019, 0.806475865866949, NA, 0.792992515529661,

0.633927820899974, NA, 0.792992515529661, 0.751416088683921,

NA, NA, -0.136965855073157, NA, NA, 0.165514438477573, NA, NA,

0.0295588022415444, NA, NA, 1.26129787094521, NA, 1.09527338740260,

1.41585316336144, NA, 0.89199803930511, 0.444685821261446, NA,

0.500775287912489, 0.182321556793955, NA, 0.506817602368452,

1.62136648329937, NA, NA, -0.211956361923645, NA, NA, 1.44691898293633,

NA, NA, -0.136965855073157, NA, NA, 3.89995042419388, NA, 3.73766961828337,

0.683096844706444, NA, 0.78845736036427, 4.05178494780330, NA,

3.27714473299218, 0.62057648772511, NA, 0.647103242058538, -0.98886142470899,

NA, NA, 1.97408102602201, NA, NA, -0.574475650842447, NA, NA,

1.60542989103656, NA, NA, -0.644357016390513, NA, -0.78965808094079,

0.683096844706444, NA, 0.741937344729377, -0.263965545834465,

NA, -0.378336440719912, 0.405465108108164, NA, 0.512823626428664,

0.357674444271816, NA, NA, -0.00904074465214907, NA, NA, 0.46373401623214,

NA, NA, -0.415515443961666, NA, NA, -0.0534007767271153, NA,

NA, -0.805196684368568, NA, NA, 0.0861776962410524, NA, NA, -0.683196849706777,

NA, NA), .Dim = c(47, 2, 2,3))

)

)

Results

                     

node         

mean        

sd             

MC error 

2.5%               

median        

97.5%         

start 

sample 

 cor[1,1]

0.5346  

0.1315

0.001055 

0.2439 

0.5469       

0.7561       

1001

20000

  cor[1,2]    

0.2567 

0.2124 

0.001507 

-0.19

0.2702

0.628

1001

  20000  

cor[2,1]   

0.1869    

  0.1755     

0.001291

-0.1721

0.1924

   0.5121 

1001

20000 

  cor[2,2]   

   0.3579

  0.1995 

     0.001492 

  -0.07114

0.3726  

0.7011 

1001 

20000

  sigmae[1] 

0.008349  

0.001687

    6.761E-5 

  0.005601 

0.00814    

0.01218 

1001

20000 

sigmae[2] 

0.006677

0.003019  

1.499E-4 

0.001425 

    0.006513       

0.01312  

 1001

20000 

h2[1]  

   0.5536

0.4989   

0.003914 

-0.4021

0.542  

1.567  

  21001 

   20000  

  h2[2]     

-0.3356

0.5348  

    0.004365   

-1.354   

-0.35    

  0.747 

    21001 

        20000