##############################
### Model to predict within-group group relatedness in mammals
### From Dyble and Clutton-Brock "Contrasts in kinship structure in mammalian societies"
### Behavioural Ecology (2020)
### Written by M. Dyble, January 2019
###########################

    ### NB: The below values are those for Damaraland mole-rats
    Nm <- 10   ## Number of adult males
    Nf <- 10   ## Number of adult females
    k <- 3    ## Litter size
    B <- 1  ## Probability of a subordinate female reproducting
    n <- 3    ## Number of juvenile cohorts at any time
    t_fem <- 0  ## Probability of a female dominant retaining her dominance from one reproductive bout to the next
    theta <- 3 ## Numer of juveniles in the group for every adult
    adult_rel <- "F"  ## Is there relatedness between adult group mates of difference sexes? Can be "Neither", "Both", "M", or "F"
    
    
    t_male <- 1 ## Probability of a male dominant retaining his dominance from one reproductive bout to the next
    alpha <- 1  ## 
    
    
    
    ### Probability of two juveniles sharing a father (Equation 1)
    new_d <- (1+ alpha*(Nm-1))/Nm
    fath <- new_d^2 + (((1-new_d)^2)/(Nm-1))
    
    temp <- 0
    if(n>1){
      for(i in 1:(n-1)){
        temp <- temp + 2*(n-i)*  (fath*(t_male^i) +  ((1-(t_male^i))/Nm))
      }}
    prob_fath <- (temp + fath*n)/(n^2)
    
    
    ### Probability of two juveniles from same cohort sharing a mother  (Equation 2)
    sum <- 0
    if(k>1){
      for(i in 1:Nf){
        sum <- sum + (choose((Nf-1),(i-1))*(k-1)* ((1-B)^(Nf-i))* (B^(i-1))) /((k*i)-1)
      }} 
    Moth <- sum
    
    
    ### Probability of two juveniles from different cohort sharing a mother i (Equation 3)
    sum <- 0
    for(i in 1:Nf){
      sum <- sum + (choose((Nf-1),(i-1))* ((1-B)^(Nf-i))* (B^(i-1))) /i
    } 
    Moth_prime <- sum
    
    ### Overall probability of juveniles sharing a mother (equation 4)
    temp <- 0
    if(n>1){for(i in 1:(n-1)){
      temp <- temp + (2*(n-i))*(((t_fem^i)  *Moth_prime) +  ((1-(t_fem^i))/Nf))
    }}
    prob_moth <- ((Moth*n) + temp)/(n^2)
    
    

    
    ### Relatedness among juveniles through shared parentage alone
    rj <- 0.25 * (prob_moth+prob_fath)
    
    ### Relatedness among juveniles and adult through shared ancestry beyond parents (Equation 6 + 7)
    if(adult_rel=="Neither") {
      ra <- 0
      rb <-  1/(Nf+Nm)
    } else
      
      if(adult_rel=="Both") {
        for(i in 1:4){rj <- (rj/4) + ((0.25-(rj/4))*prob_moth) + (rj/4) + ((0.25-(rj/4))*prob_fath)}
        ra <- ((((Nf^2)-Nf)+((Nm^2)-Nm))*rj)  / (((Nm+Nf)^2)-(Nm+Nf))   
        rb <- (1+  ((rj*(Nf-1))/2) + ((rj*(Nm-1))/2)) / (Nf+Nm) 
      } else
        
        if(adult_rel=="F"){
          for(i in 1:4){rj <- (rj/4) + ((0.25-(rj/4))*prob_moth) + (0.25*prob_fath)}
          ra <-  (((Nf^2)-Nf)*rj)  / (((Nm+Nf)^2)-(Nm+Nf))
          rb <- (1+  ((rj*(Nf-1))/2)) / (Nf+Nm)
        } else
          
          if(adult_rel=="M"){
            for(i in 1:4){rj <- (rj/4) + ((0.25-(rj/4))*prob_fath) + (0.25*prob_moth)}
            ra <-   (((Nm^2)-Nm)*rj)  / (((Nm+Nf)^2)-(Nm+Nf)) 
            rb <- (1+  ((rj*(Nm-1))/2)) / (Nf+Nm) 
          } 
    
    ### Relatedness for the whole group (equation 9)
    
    Na <- Nm+Nf
    Nj <- theta*Na
    rg <-  ((ra*Na*(Na-1))+(2*rb*Na*Nj)+(rj*Nj*((Nj)-1)))/((Na+Nj)*(Na+Nj-1))

    
    ######## 
    ## Mean relatedness in the group?
    ########
    rg