Nothing
simul.replace.uni <-
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")
#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
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
others <- data[which(data$year == generation & data$sex == 1), ]$id
others <- others[-which(others == known[i, 2])]
sire[i] <- others[sample(1:length(others), 1)]
}
#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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.