Nothing
compare.clines <-
function(cline.data1=NULL, cline.data2=NULL, sig.test=FALSE,
n.reps=1000){
## verify that input data objects were supplied
if (is.null(cline.data1)==TRUE | is.null(cline.data2)==TRUE) stop("cline.data was not supplied")
if (is.list(cline.data1)==FALSE | is.list(cline.data2)==FALSE) stop("cline.data is not a list object")
## make sure both data sets have the same number of loci
if (dim(cline.data1$Loci.data)[1] != dim(cline.data2$Loci.data)[1])
stop("the number of markers differs between data sets")
## determine if the range of values for hi.index1 > hi.index2
if (min(cline.data1$hybrid.index) > min(cline.data2$hybrid.index))
warning("regression model predictions are being extended beyond the range of the observed values used to fit the model")
if (max(cline.data1$hybrid.index) < max(cline.data2$hybrid.index))
warning("regression model predictions are being extended beyond the range of the observed values used to fit the model")
## compute likelihood ratios without significance testing
if (sig.test==TRUE) cat("compare.clines permutation test working\n this may take a few minutes\n")
if (sig.test==FALSE){
## set up data objects
n.loci<-dim(cline.data2$Loci.data)[1]
n.ind<-dim(cline.data2$Count.matrix)[2]
## compute fitted values for M1D2
M1D2.fitted.AA<-array(dim=c(n.loci,n.ind))
M1D2.fitted.Aa<-array(dim=c(n.loci,n.ind))
M1D2.fitted.aa<-array(dim=c(n.loci,n.ind))
for (i in 1:n.loci){
reg.out<-multinom(cline.data1$Count.matrix[i,]~cline.data1$hybrid.index, trace=FALSE)
## for dominant data
if (cline.data1$Loci.data[i,2]=="D" | cline.data1$Loci.data[i,2]=="d" | cline.data1$Loci.data[i,2]=="H" | cline.data1$Loci.data[i,2]=="h"){
if (length(coef(reg.out))<2) warning("warning, invariant locus included")
else{
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Hx<-exp(AA.int+AA.slope*cline.data2$hybrid.index)
M1D2.fitted.AA[i,]<-exp(AA.int+AA.slope*cline.data2$hybrid.index)/(1+Hx)
M1D2.fitted.aa[i,]<-1-M1D2.fitted.AA[i,]
}
}
## for co-dominant data
## this now uses the fit.c.clines function
else if (cline.data1$Loci.data[i,2]=="C" | cline.data1$Loci.data[i,2]=="c"){
if (length(coef(reg.out))>2){
AA.slope<-coef(reg.out)[2,2]
AA.int<-coef(reg.out)[2,1]
Aa.slope<-coef(reg.out)[1,2]
Aa.int<-coef(reg.out)[1,1]
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int) + exp(Aa.slope*cline.data2$hybrid.index+Aa.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-exp(Aa.slope*cline.data2$hybrid.index+Aa.int)/(1+Hx)
M1D2.fitted.aa[i,]<-1-(M1D2.fitted.AA[i,] + M1D2.fitted.Aa[i,])
}
else if (length(coef(reg.out))==2){
if(sum(cline.data1$Count.matrix[i,]==2,na.rm=TRUE)>=1 & sum(cline.data1$Count.matrix[i,]==1,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-1-(exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx))
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
else if(sum(cline.data1$Count.matrix[i,]==2,na.rm=TRUE)>=1 & sum(cline.data1$Count.matrix[i,]==0,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-1-(exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx))
}
else if(sum(cline.data1$Count.matrix[i,]==1,na.rm=TRUE)>=1 & sum(cline.data1$Count.matrix[i,]==0,na.rm=TRUE)>=1){
AA.slope<-NA
AA.int<-NA
Aa.slope<-coef(reg.out)[2]
Aa.int<-coef(reg.out)[1]
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int)
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-exp(Aa.slope*cline.data2$hybrid.index+Aa.int)/(1+Hx)
M1D2.fitted.aa[i,]<-1-(exp(Aa.slope*cline.data2$hybrid.index+Aa.int)/(1+Hx))
}
}
else if (length(coef(reg.out))<2){
if(sum(cline.data1$Count.matrix[i,]==2,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(1,n.ind)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(cline.data1$Count.matrix[i,]==1,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-rep(1,n.ind)
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(cline.data1$Count.matrix[i,]==0,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-rep(1,n.ind)
}
}
}
}
## compute L(M1|D2)
ln.likelihood.M1<-numeric(n.loci)
ln.likelihood.M2<-numeric(n.loci)
ln.likelihood.ratio<-numeric(n.loci)
for (i in 1:n.loci){
probM1D2<-rep(NA,n.ind)
probM2D2<-rep(NA,n.ind)
## for co-dominant markers
if (cline.data2$Loci.data[i,2]=="C" | cline.data2$Loci.data[i,2]=="c"){
for (z in 1:n.ind){
if (is.na(cline.data2$Count.matrix[i,z])==FALSE){
if (cline.data2$Count.matrix[i,z]==2){
probM1D2[z]<-M1D2.fitted.AA[i,z]
probM2D2[z]<-cline.data2$Fitted.AA[i,z]
}
else if (cline.data2$Count.matrix[i,z]==1){
probM1D2[z]<-M1D2.fitted.Aa[i,z]
probM2D2[z]<-cline.data2$Fitted.Aa[i,z]
}
else if (cline.data2$Count.matrix[i,z]==0){
probM1D2[z]<-M1D2.fitted.aa[i,z]
probM2D2[z]<-cline.data2$Fitted.aa[i,z]
}
}
else if (is.na(cline.data2$Count.matrix[i,z])==TRUE){
probM1D2[z]<-0.33
probM2D2[z]<-0.33
}
}
}
else if (cline.data2$Loci.data[i,2]=="D" | cline.data2$Loci.data[i,2]=="d" | cline.data2$Loci.data[i,2]=="H" | cline.data2$Loci.data[i,2]=="h"){
for (z in 1:n.ind){
if (is.na(cline.data2$Count.matrix[i,z])==FALSE){
if (cline.data2$Count.matrix[i,z]==2){
probM1D2[z]<-M1D2.fitted.AA[i,z]
probM2D2[z]<-cline.data2$Fitted.AA[i,z]
}
else if (cline.data2$Count.matrix[i,z]==0){
probM1D2[z]<-M1D2.fitted.aa[i,z]
probM2D2[z]<-cline.data2$Fitted.aa[i,z]
}
}
else if (is.na(cline.data2$Count.matrix[i,z])==TRUE){
probM1D2[z]<-0.5
probM2D2[z]<-0.5
}
}
}
ln.likelihood.M1[i]<-sum(log(probM1D2),na.rm=TRUE)
ln.likelihood.M2[i]<-sum(log(probM2D2),na.rm=TRUE)
ln.likelihood.ratio[i]<-ln.likelihood.M2[i]-ln.likelihood.M1[i]
}
likelihood.matrix<-array(ln.likelihood.ratio,dim=c(n.loci,1))
rownames(likelihood.matrix)<-cline.data2$Loci.data[,1]
colnames(likelihood.matrix)<-"log.likelihood.ratio"
return (likelihood.matrix)
}
## conduct permutation significance test, this is a full permutation, minimum and maximum values are not constrained
else if (sig.test==TRUE){
## set up data objects
n.loci<-dim(cline.data2$Loci.data)[1]
n.ind<-dim(cline.data2$Count.matrix)[2]
p.value<-numeric(n.loci)
## compute fitted values for M1D2
M1D2.fitted.AA<-array(dim=c(n.loci,n.ind))
M1D2.fitted.Aa<-array(dim=c(n.loci,n.ind))
M1D2.fitted.aa<-array(dim=c(n.loci,n.ind))
M2D2.fitted.AA<-array(dim=c(n.loci,n.ind))
M2D2.fitted.Aa<-array(dim=c(n.loci,n.ind))
M2D2.fitted.aa<-array(dim=c(n.loci,n.ind))
for (i in 1:n.loci){
reg.out<-multinom(cline.data1$Count.matrix[i,]~cline.data1$hybrid.index, trace=FALSE)
## for dominant data
if (cline.data1$Loci.data[i,2]=="D" | cline.data1$Loci.data[i,2]=="d" | cline.data1$Loci.data[i,2]=="H" | cline.data1$Loci.data[i,2]=="h"){
if (length(coef(reg.out))<2) warning("warning, invariant locus included")
else{
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Hx<-exp(AA.int[i]+AA.slope*cline.data2$hybrid.index)
M1D2.fitted.AA[i,]<-exp(AA.int+AA.slope*cline.data2$hybrid.index)/(1+Hx)
M1D2.fitted.aa[i,]<-1-M1D2.fitted.AA[i,]
}
}
## for co-dominant data
## this now uses the fit.c.clines function
else if (cline.data1$Loci.data[i,2]=="C" | cline.data1$Loci.data[i,2]=="c"){
if (length(coef(reg.out))>2){
AA.slope<-coef(reg.out)[2,2]
AA.int<-coef(reg.out)[2,1]
Aa.slope<-coef(reg.out)[1,2]
Aa.int<-coef(reg.out)[1,1]
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int) + exp(Aa.slope*cline.data2$hybrid.index+Aa.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-exp(Aa.slope*cline.data2$hybrid.index+Aa.int)/(1+Hx)
M1D2.fitted.aa[i,]<-1-(M1D2.fitted.AA[i,] + M1D2.fitted.Aa[i,])
}
else if (length(coef(reg.out))==2){
if(sum(cline.data1$Count.matrix[i,]==2,na.rm=TRUE)>=1 & sum(cline.data1$Count.matrix[i,]==1,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-1-(exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx))
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
else if(sum(cline.data1$Count.matrix[i,]==2,na.rm=TRUE)>=1 & sum(cline.data1$Count.matrix[i,]==0,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-1-(exp(AA.slope*cline.data2$hybrid.index+AA.int)/(1+Hx))
}
else if(sum(cline.data1$Count.matrix[i,]==1,na.rm=TRUE)>=1 & sum(cline.data1$Count.matrix[i,]==0,na.rm=TRUE)>=1){
AA.slope<-NA
AA.int<-NA
Aa.slope<-coef(reg.out)[2]
Aa.int<-coef(reg.out)[1]
Hx<-exp(AA.slope*cline.data2$hybrid.index+AA.int)
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-exp(Aa.slope*cline.data2$hybrid.index+Aa.int)/(1+Hx)
M1D2.fitted.aa[i,]<-1-(exp(Aa.slope*cline.data2$hybrid.index+Aa.int)/(1+Hx))
}
}
else if (length(coef(reg.out))<2){
if(sum(cline.data1$Count.matrix[i,]==2,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(1,n.ind)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(cline.data1$Count.matrix[i,]==1,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-rep(1,n.ind)
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(cline.data1$Count.matrix[i,]==0,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-rep(1,n.ind)
}
}
}
}
## compute L(M1|D2)
ln.likelihood.M1<-numeric(n.loci)
ln.likelihood.M2<-numeric(n.loci)
ln.likelihood.ratio<-numeric(n.loci)
for (i in 1:n.loci){
probM1D2<-rep(NA,n.ind)
probM2D2<-rep(NA,n.ind)
## for co-dominant markers
if (cline.data2$Loci.data[i,2]=="C" | cline.data2$Loci.data[i,2]=="c"){
for (z in 1:n.ind){
if (is.na(cline.data2$Count.matrix[i,z])==FALSE){
if (cline.data2$Count.matrix[i,z]==2){
probM1D2[z]<-M1D2.fitted.AA[i,z]
probM2D2[z]<-cline.data2$Fitted.AA[i,z]
}
else if (cline.data2$Count.matrix[i,z]==1){
probM1D2[z]<-M1D2.fitted.Aa[i,z]
probM2D2[z]<-cline.data2$Fitted.Aa[i,z]
}
else if (cline.data2$Count.matrix[i,z]==0){
probM1D2[z]<-M1D2.fitted.aa[i,z]
probM2D2[z]<-cline.data2$Fitted.aa[i,z]
}
}
else if (is.na(cline.data2$Count.matrix[i,z])==TRUE){
probM1D2[z]<-0.33
probM2D2[z]<-0.33
}
}
}
else if (cline.data2$Loci.data[i,2]=="D" | cline.data2$Loci.data[i,2]=="d" | cline.data2$Loci.data[i,2]=="H" | cline.data2$Loci.data[i,2]=="h"){
for (z in 1:n.ind){
if (is.na(cline.data2$Count.matrix[i,z])==FALSE){
if (cline.data2$Count.matrix[i,z]==2){
probM1D2[z]<-M1D2.fitted.AA[i,z]
probM2D2[z]<-cline.data2$Fitted.AA[i,z]
}
else if (cline.data2$Count.matrix[i,z]==0){
probM1D2[z]<-M1D2.fitted.aa[i,z]
probM2D2[z]<-cline.data2$Fitted.aa[i,z]
}
}
else if (is.na(cline.data2$Count.matrix[i,z])==TRUE){
probM1D2[z]<-0.5
probM2D2[z]<-0.5
}
}
}
ln.likelihood.M1[i]<-sum(log(probM1D2),na.rm=TRUE)
ln.likelihood.M2[i]<-sum(log(probM2D2),na.rm=TRUE)
ln.likelihood.ratio[i]<-ln.likelihood.M2[i]-ln.likelihood.M1[i]
## Permutations for neutral data set
n.total<-dim(cline.data1$Count.matrix)[2] + dim(cline.data2$Count.matrix)[2]
allCount<-cbind(cline.data1$Count.matrix,cline.data2$Count.matrix)
allHI<-c(cline.data1$hybrid.index,cline.data2$hybrid.index)
perm.ln.likelihood.ratio<-numeric(n.reps)
for(k in 1:n.reps){
## sample individual indexes for permutation
permData<-sample(1:n.total,n.total,replace=FALSE)
sam1<-permData[1:dim(cline.data1$Count.matrix)[2]]
sam2<-permData[(1 + dim(cline.data1$Count.matrix)[2]):n.total]
## Regression for first sample
reg.out<-multinom(allCount[i,sam1]~allHI[sam1], trace=FALSE)
## for dominant data
if (cline.data1$Loci.data[i,2]=="D" | cline.data1$Loci.data[i,2]=="d" | cline.data1$Loci.data[i,2]=="H" | cline.data1$Loci.data[i,2]=="h"){
if (length(coef(reg.out))<2) warning("warning, invariant locus included")
else{
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Hx<-exp(AA.int[i]+AA.slope*allHI[sam2])
M1D2.fitted.AA[i,]<-exp(AA.int+AA.slope*allHI[sam2])/(1+Hx)
M1D2.fitted.aa[i,]<-1-M1D2.fitted.AA[i,]
}
}
## for co-dominant data
else if (cline.data1$Loci.data[i,2]=="C" | cline.data1$Loci.data[i,2]=="c"){
if (length(coef(reg.out))>2){
AA.slope<-coef(reg.out)[2,2]
AA.int<-coef(reg.out)[2,1]
Aa.slope<-coef(reg.out)[1,2]
Aa.int<-coef(reg.out)[1,1]
Hx<-exp(AA.slope*allHI[sam2]+AA.int) + exp(Aa.slope*allHI[sam2]+Aa.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-exp(Aa.slope*allHI[sam2]+Aa.int)/(1+Hx)
M1D2.fitted.aa[i,]<-1-(M1D2.fitted.AA[i,] + M1D2.fitted.Aa[i,])
}
else if (length(coef(reg.out))==2){
if(sum(allCount[i,sam1]==2,na.rm=TRUE)>=1 & sum(allCount[i,sam1]==1,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*allHI[sam2]+AA.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-1-(exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx))
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
else if(sum(allCount[i,sam1]==2,na.rm=TRUE)>=1 & sum(allCount[i,sam1]==0,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*allHI[sam2]+AA.int)
M1D2.fitted.AA[i,]<-exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-1-(exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx))
}
else if(sum(allCount[i,sam1]==1,na.rm=TRUE)>=1 & sum(allCount[i,sam1]==0,na.rm=TRUE)>=1){
AA.slope<-NA
AA.int<-NA
Aa.slope<-coef(reg.out)[2]
Aa.int<-coef(reg.out)[1]
Hx<-exp(AA.slope*allHI[sam2]+AA.int)
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-exp(Aa.slope*allHI[sam2]+Aa.int)/(1+Hx)
M1D2.fitted.aa[i,]<-1-(exp(Aa.slope*allHI[sam2]+Aa.int)/(1+Hx))
}
}
else if (length(coef(reg.out))<2){
if(sum(allCount[i,sam1]==2,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(1,n.ind)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(allCount[i,sam1]==1,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-rep(1,n.ind)
M1D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(allCount[i,sam1]==0,na.rm=TRUE)>=1){
M1D2.fitted.AA[i,]<-rep(0,n.ind)
M1D2.fitted.Aa[i,]<-rep(0,n.ind)
M1D2.fitted.aa[i,]<-rep(1,n.ind)
}
}
## compute L(M1|D2)
probM1D2<-rep(NA,n.ind)
## for dominant markers
if (cline.data1$Loci.data[i,2]=="D" | cline.data1$Loci.data[i,2]=="d" | cline.data1$Loci.data[i,2]=="H" | cline.data1$Loci.data[i,2]=="h"){
for (z in 1:n.ind){
zSam2<-sam2[z]
if (is.na(allCount[i,zSam2]==FALSE)){
if (allCount[i,zSam2]==2) probM1D2[z]<-M1D2.fitted.AA[i,z]
else if (allCount[i,zSam2]==0) probM1D2[z]<-M1D2.fitted.aa[i,z]
}
else if (is.na(allCount[i,zSam2])==TRUE) probM1D2[z]<-0.5
}
}
## for co-dominant markers
else if (cline.data1$Loci.data[i,2]=="C" | cline.data1$Loci.data[i,2]=="c"){
for (z in 1:n.ind){
zSam2<-sam2[z]
if (is.na(allCount[i,zSam2])==FALSE){
if (allCount[i,zSam2]==2) probM1D2[z]<-M1D2.fitted.AA[i,z]
else if (allCount[i,zSam2]==1) probM1D2[z]<-M1D2.fitted.Aa[i,z]
else if (allCount[i,zSam2]==0) probM1D2[z]<-M1D2.fitted.aa[i,z]
}
else if (is.na(allCount[i,zSam2])==TRUE) probM1D2[z]<-0.33
}
}
replnLikM1D2<-sum(log(probM1D2))
## Now calculate M2D2
## Regression for first sample
reg.out<-multinom(allCount[i,sam2]~allHI[sam2], trace=FALSE)
if (cline.data1$Loci.data[i,2]=="D" | cline.data1$Loci.data[i,2]=="d" | cline.data1$Loci.data[i,2]=="H" | cline.data1$Loci.data[i,2]=="h"){
if (length(coef(reg.out))<2) warning("warning, invariant locus included")
else{
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Hx<-exp(AA.int[i]+AA.slope*allHI[sam2])
M2D2.fitted.AA[i,]<-exp(AA.int+AA.slope*allHI[sam2])/(1+Hx)
M2D2.fitted.aa[i,]<-1-M2D2.fitted.AA[i,]
}
}
else if (cline.data1$Loci.data[i,2]=="C" | cline.data1$Loci.data[i,2]=="c"){
if (length(coef(reg.out))>2){
AA.slope<-coef(reg.out)[2,2]
AA.int<-coef(reg.out)[2,1]
Aa.slope<-coef(reg.out)[1,2]
Aa.int<-coef(reg.out)[1,1]
Hx<-exp(AA.slope*allHI[sam2]+AA.int) + exp(Aa.slope*allHI[sam2]+Aa.int)
M2D2.fitted.AA[i,]<-exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx)
M2D2.fitted.Aa[i,]<-exp(Aa.slope*allHI[sam2]+Aa.int)/(1+Hx)
M2D2.fitted.aa[i,]<-1-(M2D2.fitted.AA[i,] + M2D2.fitted.Aa[i,])
}
else if (length(coef(reg.out))==2){
if(sum(allCount[i,sam2]==2,na.rm=TRUE)>=1 & sum(allCount[i,sam2]==1,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*allHI[sam2]+AA.int)
M2D2.fitted.AA[i,]<-exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx)
M2D2.fitted.Aa[i,]<-1-(exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx))
M2D2.fitted.aa[i,]<-rep(0,n.ind)
}
else if(sum(allCount[i,sam2]==2,na.rm=TRUE)>=1 & sum(allCount[i,sam2]==0,na.rm=TRUE)>=1){
AA.slope<-coef(reg.out)[2]
AA.int<-coef(reg.out)[1]
Aa.slope<-NA
Aa.int<-NA
Hx<-exp(AA.slope*allHI[sam2]+AA.int)
M2D2.fitted.AA[i,]<-exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx)
M2D2.fitted.Aa[i,]<-rep(0,n.ind)
M2D2.fitted.aa[i,]<-1-(exp(AA.slope*allHI[sam2]+AA.int)/(1+Hx))
}
else if(sum(allCount[i,sam2]==1,na.rm=TRUE)>=1 & sum(allCount[i,sam2]==0,na.rm=TRUE)>=1){
AA.slope<-NA
AA.int<-NA
Aa.slope<-coef(reg.out)[2]
Aa.int<-coef(reg.out)[1]
Hx<-exp(AA.slope*allHI[sam2]+AA.int)
M2D2.fitted.AA[i,]<-rep(0,n.ind)
M2D2.fitted.Aa[i,]<-exp(Aa.slope*allHI[sam2]+Aa.int)/(1+Hx)
M2D2.fitted.aa[i,]<-1-(exp(Aa.slope*allHI[sam2]+Aa.int)/(1+Hx))
}
}
else if (length(coef(reg.out))<2){
if(sum(allCount[i,sam2]==2,na.rm=TRUE)>=1){
M2D2.fitted.AA[i,]<-rep(1,n.ind)
M2D2.fitted.Aa[i,]<-rep(0,n.ind)
M2D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(allCount[i,sam2]==1,na.rm=TRUE)>=1){
M2D2.fitted.AA[i,]<-rep(0,n.ind)
M2D2.fitted.Aa[i,]<-rep(1,n.ind)
M2D2.fitted.aa[i,]<-rep(0,n.ind)
}
if(sum(allCount[i,sam2]==0,na.rm=TRUE)>=1){
M2D2.fitted.AA[i,]<-rep(0,n.ind)
M2D2.fitted.Aa[i,]<-rep(0,n.ind)
M2D2.fitted.aa[i,]<-rep(1,n.ind)
}
}
}
## compute L(M2|D2)
probM2D2<-rep(NA,n.ind)
## for dominant markers
if (cline.data1$Loci.data[i,2]=="D" | cline.data1$Loci.data[i,2]=="d" | cline.data1$Loci.data[i,2]=="H" | cline.data1$Loci.data[i,2]=="h"){
for (z in 1:n.ind){
zSam2<-sam2[z]
if (is.na(allCount[i,zSam2]==FALSE)){
if (allCount[i,zSam2]==2) probM2D2[z]<-M2D2.fitted.AA[i,z]
else if (allCount[i,zSam2]==0) probM2D2[z]<-M2D2.fitted.aa[i,z]
}
else if (is.na(allCount[i,zSam2])==TRUE) probM2D2[z]<-0.5
}
}
## for co-dominant markers
else if (cline.data1$Loci.data[i,2]=="C" | cline.data1$Loci.data[i,2]=="c"){
for (z in 1:n.ind){
zSam2<-sam2[z]
if (is.na(allCount[i,zSam2])==FALSE){
if (allCount[i,zSam2]==2) probM2D2[z]<-M2D2.fitted.AA[i,z]
else if (allCount[i,zSam2]==1) probM2D2[z]<-M2D2.fitted.Aa[i,z]
else if (allCount[i,zSam2]==0) probM2D2[z]<-M2D2.fitted.aa[i,z]
}
else if (is.na(allCount[i,zSam2])==TRUE) probM2D2[z]<-0.33
}
}
replnLikM2D2<-sum(log(probM2D2))
perm.ln.likelihood.ratio[k]<-replnLikM2D2 - replnLikM1D2
}
## This ends the rep loop
}
p.value[i]<-sum(perm.ln.likelihood.ratio > ln.likelihood.ratio[i])/n.reps
cat("permutation test for locus: ", cline.data1$Loci.data[i,1],"\n")
## This ends the locus loop
}
## This ends the locus loop
likelihood.matrix<-cbind(ln.likelihood.ratio,p.value)
rownames(likelihood.matrix)<-cline.data2$Loci.data[,1]
colnames(likelihood.matrix)<-c("log.likelihood.ratio","PValues")
return (likelihood.matrix)
}
}
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.