R/simul.replace.similar.R

In PSIMEX: SIMEX Algorithm on Pedigree Structures

simul.replace.similar <-
function(pedigree, lambda, B,   data, model) {
  
  
  
  #where to generate more errors
  known <- pedigree[which(!is.na(pedigree$dam)), ]
  #already missing
  unknown <- pedigree[which(is.na(pedigree$dam)), ]
 
  names(known) <- c("animal", "sire", "dam")
  names(unknown) <- c("animal", "sire", "dam")
  resp <- all.vars(model$call)[1]
  

 
  #setting initial values
  inb <- matrix(data  =  NA, nrow  =  length(lambda), ncol  =  B)
  se_inb <- matrix(data  =  NA, nrow  =  length(lambda), ncol  =  B)
  mean_inb <- matrix(data  =  NA, nrow  =  length(lambda), ncol = B)
  median_inb <- matrix(data  =  NA, nrow  =  length(lambda), ncol = B)
  var_inb <- matrix(data  =  NA, nrow  =  length(lambda), ncol = B)
  pval <- matrix(data  =  NA, nrow  =  length(lambda), ncol = B)
  
  for (k in 1:length(lambda)){
    for (j in 1:B) {
      animal <- known$animal
      dam <- known$dam
      sire <- known$sire  
      
      #initializing inbreeding coefficient
      data$f_inb<- c()
      
      #choosing lambda over the total number of individuals 
      #and setting their parents as random
      m <- sample(1:length(known$animal), lambda[k]*length(known$animal))
      
      for (i in m) {
        
        generation <- data[which(data$id == known$sire[i]), ]$year
        dad <- pedigree[which(pedigree$animal == known$animal[i]), 2]
        others <- data[which(data$year == generation & data$sex == 1),  ]$id
        others <- others[-which(others == known[i, 2])]
        #compute the difference among the dad's trait and all the others'
        trait_dad <- data[data$id == dad, which(names(data) == resp)]
        trait_others <- data[which(is.element(data$id, others)), which(names(data) == resp)]
        diff <- trait_dad - trait_others
        #choose the most similar one
        similar <- which(abs(diff) == min(abs(diff), na.rm = TRUE))
        sire[i] <- others[similar]
        
      }
      
      
      #recalculate pedigree
      ped <- data.frame(animal = animal, dam = dam, sire = sire)
      ped <- rbind(ped, unknown[c('animal', 'dam', 'sire')])
      ord <- orderPed(ped)
      ped <- ped[order(ord),]
      
      #re calculate inbreeding coefficient
      f_inb <- data.frame(ped$animal, calcInbreeding(ped))
      names(f_inb) <- c("id", "f_inb")
      data <- merge(data, f_inb, by = "id")
      data$f_inb <- as.numeric(data$f_inb)
      
      if(sum(data$f_inb)>0){
        
        #model 
        #inbreeding depression is the regression coefficient of f_inb
        type <- as.character(model$call[1])
        if (type == 'lm'){
          model1 <- lm(formula = as.character(model$call[2]), data = data)
        }
        
        if (type == 'glm'){
          model1 <- glm(formula = as.character(model$call[2]), family = model$call$family,  data = data)
        }
        
        
        #inbreeding depression
        inb[k, j] <- summary(model1)$coefficients["f_inb", 1]
        se_inb[k,j] <- summary(model1)$coefficients["f_inb", 2]
        pval [k, j] <- summary(model1)$coefficients["f_inb", 4]
        
        #info on inbreeding coefficient
        mean_inb[k, j] <- mean(data$f_inb)
        median_inb[k, j] <- median(data$f_inb)
        var_inb[k, j] <- var(data$f_inb)
        
      }
      
      
      #this must be checked
      if(sum(data$f_inb) == 0) {
        inb[k, j] <- NA
        se_inb[k,j] <- NA
        mean_inb[k, j] <- NA
        median_inb[k, j] <- NA
        var_inb[k, j] <- NA
        pval [k, j] <- NA
        
      }
      
    }
    
  }
  
  #store results
  res = list(inb, se_inb, pval, 
            mean_inb, median_inb, var_inb)
  names(res) <- c('inb_dep', 'se_inb_dep', 'pval_inb_dep', 'mean_inb', 'median_inb', 'var_inb')
  return(res)
}