simX {BookEKM} | R Documentation |
Data for pairwise relationships are simulated. LR is calculated and also summary statistics (mean, SD). For X, E(LR|numerator)=1-k1^2+k1^2*(L+1)/2
simX(x1, id1, id2, p, nsim, X = TRUE)
x1 |
Linkdat (pedigree) |
id1 |
id of female |
id2 |
id of female |
p |
allele frequency |
nsim |
number of simulations |
X |
|
mean, SD and LR-s
Thore.Egeland@nmbu.no
## Not run: library(paramlink) ### Example 1 # Table 1 of mdvX.pdf note with L=10 set.seed(17) XMean <- AutoMean <- NULL L <- 10 p <- runif(L); p <- p/sum(p) nsim <- 1000 x1 <- nuclearPed(1) id1 <- 2; id2 <- 3 XMean <- c(XMean, simX(x1,id1,id2,p,nsim=nsim)$meanLR) AutoMean <- c(AutoMean, simX(x1,id1,id2,p,nsim=nsim, X = FALSE)$meanLR) #HSI (Half Sib Index) x1 <- nuclearPed(1) x1 <- addOffspring(x1,mother=2,noff=1,sex=2) id1 <- 5; id2 <- 3 XMean <- c(XMean, simX(x1,id1,id2,p,nsim=nsim)$meanLR) AutoMean <- c(AutoMean, simX(x1,id1,id2,p,nsim=nsim, X = FALSE)$meanLR) # maternal FC x1<-nuclearPed(2,sex=2) x1 <- addOffspring(x1,mother=3,noff=1,sex=2) x1 <- addOffspring(x1,mother=4,noff=1,sex=1) id1 <- 6; id2 <- 8 XMean <- c(XMean, simX(x1,id1,id2,p,nsim=nsim)$meanLR) AutoMean <- c(AutoMean, simX(x1,id1,id2,p,nsim=nsim, X = FALSE)$meanLR) data.frame(meanX = XMean, meanAuto = AutoMean) ### Example 2 #HSI (Half Sister Index, two markers) library(paramlink) x1 <- nuclearPed(1,sex=2) x1 <- addOffspring(x1,mother=2,noff=1,sex=2) id1 <- 5; id2 <- 3 g <- list(c(1,1),c(2,2),c(1,2)) p1 <- c(0.1, 0.9) p2 <- c(0.2, 0.8) ELR <- 0 emptySNP1 <- marker(x1, alleles=1:2, afreq=c(0.1, 0.9)) emptySNP2 <- marker(x1, alleles=1:2, afreq=c(0.2, 0.8)) t2 <- twoMarkerDistribution(x1, id=id1, emptySNP1, emptySNP2, theta=0.1) for (i in g) for (j in g){ factor <- 0.5+ (i[1]==i[2])&(j[1]==j[2]) SNP1 <- marker(x1, id2, g[[i]], alleles=1:2, afreq=c(0.1, 0.9)) SNP2 <- marker(x1, id2, g[[j]], alleles=1:2, afreq=c(0.2, 0.8)) t1 <- twoMarkerDistribution(x1, id=id1, SNP1, SNP2, theta=0.1) ELR <- ELR+factor*sum(t1^2/t2*prod(p1[g[[i]]])*prod(p2[g[[j]]])) } ### Example 3 data(NorwegianFrequencies) dd <- NorwegianFrequencies$SE33 p <- as.double(dd); p <- p/sum(p) LRx <- simX(x1,id1,id2,p,nsim=nsim) LRauto <-simX(x1,id1,id2,p,nsim=nsim,X=FALSE) #FC x1<-nuclearPed(2,sex=2) x1 <- addOffspring(x1,mother=3,noff=1,sex=2) x1 <- addOffspring(x1,mother=4,noff=1,sex=1) id1 <- 6 id2 <- 8 p <- runif(10); p <- p/sum(p) foo <- simX(x1,id1,id2,p,nsim=nsim,X=FALSE) ## End(Not run)