##sample metadata The metadata has some animal information, here we only present a few variables

meta_data
## # A tibble: 794 x 7
##    Animal Sex   Rep   Finisher_Wt Finisher_Pen Finisher_Post_F~
##     <dbl> <chr> <fct>       <dbl> <fct>                   <dbl>
##  1    107 g     2            52   11_2                     3.37
##  2     85 g     2            66   11_2                     4.49
##  3    113 g     2            61.4 11_2                     4.49
##  4    111 g     2            73.8 11_2                     4.96
##  5     79 g     2            59.4 11_2                     4.29
##  6    103 g     2            69.6 11_2                     4.33
##  7     99 g     2            53.2 11_2                     4.42
##  8     89 g     2            71.2 11_2                     3.95
##  9     55 g     2            67   11_2                     4.01
## 10     25 g     2            54   11_2                     4.78
## # ... with 784 more rows, and 1 more variable: total_attack <dbl>

Example of social interactions

The traits in the Zc matrix are time spent attacking another animal. Let’s select a few pens and look at the interaction matrices. The columns correspond to the ID of the animal delivering the attack, the rows are the animals receiving the attack. The cells contain the total time (in seconds) that each animal spent attacking each group mate over a 6 hour period of observation post-mixing.

p1<-meta_data$Animal[meta_data$Finisher_Pen=="11_2"]
p2<-meta_data$Animal[meta_data$Finisher_Pen=="10_9"]
p3<-meta_data$Animal[meta_data$Finisher_Pen=="6_3"]

Zc[as.character(p1),as.character(p1)]
##     107 85 113 111 79 103 99 89  55 25 71 69 83 27 53
## 107   0  1   1   0  7  11  6  6   5  1  5 11 14  0  2
## 85   39  0   0   4 49  31 23  4   0 21 53  7  7  0  0
## 113   5  0   0   8 65  22  1  0   3 16 13 62 12  4  0
## 111   2  0  12   0  9   0  1  6   1 11  0  2  1  8  4
## 79    1  5   6  30  0  20  5  2 109  1  1 42  6  1  0
## 103   2  5   4   1  9   0  2  4   3  3  2  7  2  0  3
## 99    1  0   0   0  3   1  0  5  24  0  0  6  2 38  4
## 89    2  4   5   5 31   0  5  0  61  6  0  0  0  0  1
## 55    0  0   0   0  0   3  1  0   0  4  0 28  6  6  4
## 25    0  9   0  22  8   1  4  5  69  0  1  1  0  0  2
## 71    0  0   0   0  0   2  0  0   0  3  0  1  2  0  1
## 69    0  1   5  33  0   1  4  1  23  2  0  0  1  0  0
## 83    1  1   1   0  1   1  3  1   0  0  4  1  0  0  1
## 27    1  0   0   0  0   0  4  2   4  0  0  0  1  0  6
## 53    1 12   0   4  8   0  0  0   1  0 17  0  1  0  0
Zc[as.character(p2),as.character(p2)]
##      1031 1075 1137 1097 1077 1200 1149 1021 1143 1023 1051 1119 1111
## 1031    0    7    0    0    0    5    4    0    2    1    0    0    1
## 1075    0    0    0    4    0    0    0    0    0    0    0    0    0
## 1137    0    0    0    0    0    2    3    7    8    1    4   11    0
## 1097   18   11   23    0   10   29    2    0    1    0    0    0    3
## 1077   40    1    0    0    0    2    0    3    0    5    7    0    3
## 1200    0    0    0    0    0    0    0    0    0    0    1    0    5
## 1149   23    0    1   17    0    6    0    0    0    3    0    0    2
## 1021   14    0    1    3    0    0    0    0    4    2    2    0    1
## 1143   23    2  196    0    3   15    7   41    0   32   17   31    9
## 1023    7    0    1    1    0    1    0    7    0    0    7    0    0
## 1051    0    0    5    0    3    0    0    2    0    4    0   18    6
## 1119    0    1    0    1    0    1    0    0    7    6    1    0    0
## 1111    0    0    0    0    0    0    0    0    5    0    1    0    0
Zc[as.character(p3),as.character(p3)]
##     139 137 183 133 195 239 243 191 245 251 131 215 277 129 213
## 139   0   3   0   0   0   0   0   0   0   0   1   0  18   0   0
## 137   2   0  12  14   1   0   2   0   1   2  24  46   6   5   1
## 183   1   6   0   1   0   0   2   3   0   0   5   4   1   0   0
## 133   0   0   3   0  12 499  20   6   1   0   2   0   2   1   2
## 195   4  32  26   3   0 241  17  66  27  34   0  60   5  26  29
## 239   3   0   1  10   6   0   0   8   0   0   0   0   3   0   5
## 243   0   2   3   2   0   0   0   2   3   4  11   6   4   4   2
## 191   6   4   0  28   7 238   0   0  17   3   3   3   6   7   7
## 245   0   0   7   0   0   0   0   0   0   1   0   0  37  10   0
## 251   0   1   1   0   1   0   1   0   2   0   1   2   0   2   1
## 131   7   3   7   1  28   2   3   5   8   2   0   2   4  11   3
## 215  10   0   1   4   4   1   0   0   0   1   3   0  13   6   7
## 277  30   0   3   4   1   0   1   1  14  14  16   3   0   1  23
## 129   0   0   0   0   0   0   0   0   0   1  10   1   0   0   1
## 213   1   0   6   2   1 345   0   0  19  10   0   3   1   0   0

Notice that in the last example there is animal (ID 239) that spend a large ammount of time attacking a few others. This animal alone spent a total of over 1300 seconds attacking other animals.

summary of the total time spent in attacks

The summaries below are computed with and without the most extreme animal in the set.

summary(meta_data$total_attack)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   36.00   61.00   75.58   97.00 1326.00
qqnorm(meta_data$total_attack)

ggplot(meta_data,aes(x=Finisher_Pen,y=total_attack))+geom_boxplot()

summary(meta_data$total_attack[-which.max(meta_data$total_attack)])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0      36      61      74      97     331
qqnorm(meta_data$total_attack[-which.max(meta_data$total_attack)])

filter(meta_data,total_attack!=max(total_attack))%>%ggplot(aes(x=Finisher_Pen,y=total_attack))+geom_boxplot()

Marginal response vs dyadic response

The total attack time is a marginal response for an individual that is the sum of all attacks to their group mates. The dyadic phenotype is the actual time spent by animal A attacking animal B.

Here we present as an example two animals with identical marginal phenotypes and very different distribution of dyadic phenotypes.

ft<-filter(meta_data,total_attack==103)
id<-as.character(ft$Animal)
pen<-as.character(ft$Finisher_Pen)
ft
## # A tibble: 2 x 7
##   Animal Sex   Rep   Finisher_Wt Finisher_Pen Finisher_Post_F~ total_attack
##    <dbl> <chr> <fct>       <dbl> <fct>                   <dbl>        <dbl>
## 1     45 g     2            64.6 3_2                      4.25          103
## 2    732 b     6            57   10_6                     3.47          103
Zc[filter(meta_data,Finisher_Pen==pen[1])%>%select(Animal)%>%unlist()%>%as.character(),id[1]]
##  49  43 101 115  51  73  67  87  45  81  47 105  93  95 109 
##   0   0   0   0   9   0  13  23   0   0   0  46   0   3   9
Zc[filter(meta_data,Finisher_Pen==pen[2])%>%select(Animal)%>%unlist()%>%as.character(),id[2]]
## 638 618 730 652 624 700 646 610 732 612 656 666 620 614 
##  13   0   1   3  23   6  21   4   0  11   5  11   3   2

Reduce all variation in Zc to 0-1?

Zi<-1*(Zc>0)
meta_data$num_attack<-colSums(Zi)

Zi[as.character(p1),as.character(p1)]
##     107 85 113 111 79 103 99 89 55 25 71 69 83 27 53
## 107   0  1   1   0  1   1  1  1  1  1  1  1  1  0  1
## 85    1  0   0   1  1   1  1  1  0  1  1  1  1  0  0
## 113   1  0   0   1  1   1  1  0  1  1  1  1  1  1  0
## 111   1  0   1   0  1   0  1  1  1  1  0  1  1  1  1
## 79    1  1   1   1  0   1  1  1  1  1  1  1  1  1  0
## 103   1  1   1   1  1   0  1  1  1  1  1  1  1  0  1
## 99    1  0   0   0  1   1  0  1  1  0  0  1  1  1  1
## 89    1  1   1   1  1   0  1  0  1  1  0  0  0  0  1
## 55    0  0   0   0  0   1  1  0  0  1  0  1  1  1  1
## 25    0  1   0   1  1   1  1  1  1  0  1  1  0  0  1
## 71    0  0   0   0  0   1  0  0  0  1  0  1  1  0  1
## 69    0  1   1   1  0   1  1  1  1  1  0  0  1  0  0
## 83    1  1   1   0  1   1  1  1  0  0  1  1  0  0  1
## 27    1  0   0   0  0   0  1  1  1  0  0  0  1  0  1
## 53    1  1   0   1  1   0  0  0  1  0  1  0  1  0  0
Zi[as.character(p2),as.character(p2)]
##      1031 1075 1137 1097 1077 1200 1149 1021 1143 1023 1051 1119 1111
## 1031    0    1    0    0    0    1    1    0    1    1    0    0    1
## 1075    0    0    0    1    0    0    0    0    0    0    0    0    0
## 1137    0    0    0    0    0    1    1    1    1    1    1    1    0
## 1097    1    1    1    0    1    1    1    0    1    0    0    0    1
## 1077    1    1    0    0    0    1    0    1    0    1    1    0    1
## 1200    0    0    0    0    0    0    0    0    0    0    1    0    1
## 1149    1    0    1    1    0    1    0    0    0    1    0    0    1
## 1021    1    0    1    1    0    0    0    0    1    1    1    0    1
## 1143    1    1    1    0    1    1    1    1    0    1    1    1    1
## 1023    1    0    1    1    0    1    0    1    0    0    1    0    0
## 1051    0    0    1    0    1    0    0    1    0    1    0    1    1
## 1119    0    1    0    1    0    1    0    0    1    1    1    0    0
## 1111    0    0    0    0    0    0    0    0    1    0    1    0    0
Zi[as.character(p3),as.character(p3)]
##     139 137 183 133 195 239 243 191 245 251 131 215 277 129 213
## 139   0   1   0   0   0   0   0   0   0   0   1   0   1   0   0
## 137   1   0   1   1   1   0   1   0   1   1   1   1   1   1   1
## 183   1   1   0   1   0   0   1   1   0   0   1   1   1   0   0
## 133   0   0   1   0   1   1   1   1   1   0   1   0   1   1   1
## 195   1   1   1   1   0   1   1   1   1   1   0   1   1   1   1
## 239   1   0   1   1   1   0   0   1   0   0   0   0   1   0   1
## 243   0   1   1   1   0   0   0   1   1   1   1   1   1   1   1
## 191   1   1   0   1   1   1   0   0   1   1   1   1   1   1   1
## 245   0   0   1   0   0   0   0   0   0   1   0   0   1   1   0
## 251   0   1   1   0   1   0   1   0   1   0   1   1   0   1   1
## 131   1   1   1   1   1   1   1   1   1   1   0   1   1   1   1
## 215   1   0   1   1   1   1   0   0   0   1   1   0   1   1   1
## 277   1   0   1   1   1   0   1   1   1   1   1   1   0   1   1
## 129   0   0   0   0   0   0   0   0   0   1   1   1   0   0   1
## 213   1   0   1   1   1   1   0   0   1   1   0   1   1   0   0
summary(meta_data$num_attack)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   5.000   7.000   7.096   9.000  13.000
qqnorm(meta_data$num_attack)

ggplot(meta_data,aes(x=Finisher_Pen,y=num_attack))+geom_boxplot()

Questions

  1. How much of the variation in marginal phenotype is explained by genetics (additive if you wish): easy answer using gblup or similar model.

  2. what are the specific genetic markers that explain the most variation in marginal phenotypes: perform a GWA

  3. What are the genetic factors that explain variation in dyadic phenotypes?

A mixed model of deliverer (ID2) and victim effect (ID1) of attacks

ID<-rownames(Zc)
ID1=rep(ID,each=length(ID))

ID2=rep(ID,times=length(ID))
y=as.vector(Zc)
length(y)
## [1] 630436
ind1<-match(ID1,meta_data$Animal)
ind2<-match(ID2,meta_data$Animal)
pen1<-meta_data$Finisher_Pen[ind1]
pen2<-meta_data$Finisher_Pen[ind2]
tmp=which((ID1==ID2)|(pen1!=pen2))
ID1=ID1[-tmp]
ID2=ID2[-tmp]
y=y[-tmp]
pen1<-pen1[-tmp]
pen2<-pen2[-tmp]
length(y)
## [1] 9958
dt=data.frame(y=y,ID1=ID1,ID2=ID2,pen1=pen1,pen2=pen2,dyad=paste(ID1,ID2,sep="_"))
dt$dyad=apply(dt,1, function(x) paste(sort(as.character(x[c(2,3)])),collapse="_"))

library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:tidyr':
## 
##     expand
lmr<-lmer(y~(1|ID1)+(1|ID2),data=dt)
lmr
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ (1 | ID1) + (1 | ID2)
##    Data: dt
## REML criterion at convergence: 81470.23
## Random effects:
##  Groups   Name        Std.Dev.
##  ID1      (Intercept)  3.496  
##  ID2      (Intercept)  5.553  
##  Residual             13.509  
## Number of obs: 9958, groups:  ID1, 794; ID2, 794
## Fixed Effects:
## (Intercept)  
##       6.055
res<-attr(VarCorr(lmr), "sc")^2
1-res/var(y)
## [1] 0.1918362
lmr<-lmer(y~(1|ID1)+(1|ID2)+(1|pen1),data=dt)
lmr
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ (1 | ID1) + (1 | ID2) + (1 | pen1)
##    Data: dt
## REML criterion at convergence: 81467.86
## Random effects:
##  Groups   Name        Std.Dev.
##  ID1      (Intercept)  3.454  
##  ID2      (Intercept)  5.475  
##  pen1     (Intercept)  1.170  
##  Residual             13.513  
## Number of obs: 9958, groups:  ID1, 794; ID2, 794; pen1, 59
## Fixed Effects:
## (Intercept)  
##       6.066
res<-attr(VarCorr(lmr), "sc")^2
1-res/var(y)
## [1] 0.1914274
lmr<-lmer(y~(1|ID1)+(1|ID2)+(1|pen1)+(1|dyad),data=dt)
lmr
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ (1 | ID1) + (1 | ID2) + (1 | pen1) + (1 | dyad)
##    Data: dt
## REML criterion at convergence: 81448.61
## Random effects:
##  Groups   Name        Std.Dev.
##  dyad     (Intercept)  3.483  
##  ID2      (Intercept)  5.469  
##  ID1      (Intercept)  3.445  
##  pen1     (Intercept)  1.142  
##  Residual             13.059  
## Number of obs: 9958, groups:  dyad, 4979; ID2, 794; ID1, 794; pen1, 59
## Fixed Effects:
## (Intercept)  
##       6.065
res<-attr(VarCorr(lmr), "sc")^2
1-res/var(y)
## [1] 0.2447491
lmr<-lmer(y~(1|ID1)+(1|ID2)+(1|dyad),data=dt)
lmr
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ (1 | ID1) + (1 | ID2) + (1 | dyad)
##    Data: dt
## REML criterion at convergence: 81450.73
## Random effects:
##  Groups   Name        Std.Dev.
##  dyad     (Intercept)  3.494  
##  ID2      (Intercept)  5.541  
##  ID1      (Intercept)  3.483  
##  Residual             13.054  
## Number of obs: 9958, groups:  dyad, 4979; ID2, 794; ID1, 794
## Fixed Effects:
## (Intercept)  
##       6.054
res<-attr(VarCorr(lmr), "sc")^2
1-res/var(y)
## [1] 0.2453903

a social relation model

Implemented following Behavior Research Methods volume 44, pages455–470(2012). These are either MLM or SEM, so we can consider these too. A limitation is that it uses a method of moments to obtain variance estimates.

Note: error and relationship variance can’t be separated because there is only one observation from dyad.

Note 2: covariances are estimated.

Actor = receives aggression Partner = delivers aggression

library(TripleR)
SRM<-RR(y~ID1*ID2|pen1,data=dt)
SRM
## Round-Robin object ('RR'), calculated by TripleR
## ------------------------------------------------
## Univariate analysis of one round robin variable in multiple groups (significance test based on Lashley & Bond, 1997, Psychological Methods)
## 
## Univariate analyses for: y 
## ---------
## Group descriptives: n =  59 ; average group size =  13.46 ; range:  10 - 15 
##                          estimate standardized    se t.value p.value
## actor variance             11.195        0.050 3.742   2.992   0.001
## partner variance           30.407        0.136 2.979  10.208   0.000
## relationship variance     181.458        0.813 4.783  37.935   0.000
## error variance                 NA           NA    NA      NA      NA
## actor-partner covariance    0.886        0.048 2.026   0.437   0.662
## relationship covariance    11.924        0.066 4.783   2.493   0.013
## Actor effect reliability: .432 
## Partner effect reliability: .674 
## NULL
1-SRM$varComp[3,2]/var(y)
## [1] 0.1964367

Observe: actor and partner variances are very close to VCE from our mixed model, relathinship variance is almost identical to residual variance in our model.

The only differences are due to covariances modeled in tripleR package