###Code for computing the optimal additional covariates for "Unbiased and Locally Efficient Estimation of Genetic Effect on Quantitative Trait ###in the Presence of Population Admixture" ################ Transmission probabilities ########################################################################################################## z2kid<-matrix(c(1,0,0,0,0,0,0,0,0, 1/4,1/4,0,1/4,1/4,0,0,0,0, 0,0,0,0,1,0,0,0,0, 1/16,1/8,1/16,1/8,1/4,1/8,1/16,1/8,1/16, 0,0,0,0,1/4,1/4,0,1/4,1/4, 0,0,0,0,0,0,0,0,1), byrow=FALSE,nrow=9) #x2kid<-t(matrix(c(2,2,2,1,2,0,1,2,1,1,1,0,0,2,0,1,0,0), nrow=9, byrow=TRUE)) #additive genetic effect x2kid<-t(matrix(c(1,1,1,1,1,0,1,1,1,1,1,0,0,1,0,1,0,0), nrow=9, byrow=TRUE)) #dominant genetic effect ################ Function to compute the optimal additional covariates to be included in the linear model ############################################### optimal<-function(p1,p2) { p<-p1 q<-1-p w<-c(p^4,4*p^3*q,2*p^2*q^2,4*p^2*q^2,4*p*q^3,q^4) wmix<-diag(c(z2kid%*%w)) miss2_mix<-x2kid%*%z2kid%*%solve(t(z2kid)%*%solve(wmix)%*%z2kid)%*%t(z2kid)%*%solve(wmix) cDD<-c(1,2,4,5) xDD<-x2kid[,cDD]; z2DD<-z2kid[cDD,1:3] wDDv<-c(p^2,2*p*q,q^2); # Pr(g*|g**) wDD<-diag(c(z2DD%*%wDDv)) #sum_{g*}(Pr(g|g*)Pr(g*|g**)=Pr(g|g**) miss1_mix_DD0<-xDD%*%z2DD%*%solve(t(z2DD)%*%solve(wDD)%*%z2DD)%*%t(z2DD)%*%solve(wDD) miss1_mix_DD<-cbind(miss1_mix_DD0[,1:2],c(0,0),miss1_mix_DD0[,3:4],c(0,0),c(0,0),c(0,0),c(0,0)) wDd<-diag(c(z2kid[,c(2,4,5)]%*%wDDv)) miss1_mix_Dd<-x2kid%*%z2kid[,c(2,4,5)]%*%solve(t(z2kid[,c(2,4,5)])%*%solve(wDd)%*%z2kid[,c(2,4,5)])%*%t(z2kid[,c(2,4,5)])%*%solve(wDd) cdd<-c(5,6,8,9) xdd<-x2kid[,cdd]; z2dd<-z2kid[cdd,c(3,5,6)] wdd<-diag(c(z2dd%*%wDDv)) miss1_mix_dd0<-xdd%*%z2dd%*%solve(t(z2dd)%*%solve(wdd)%*%z2dd)%*%t(z2dd)%*%solve(wdd) miss1_mix_dd<-cbind(c(0,0),c(0,0),c(0,0),c(0,0),miss1_mix_dd0[,1:2],c(0,0),miss1_mix_dd0[,3:4]) covamis<-array(0,c(nrep,n,2)) #matrix of the orginial covariate and new covariates to be returned for (i in 1:nrep) {for (j in 1:10) #missing 2 parents pop1 {indi<-covindi(genokid[i,(j*2-1),1:2],genokid[i,(j*2),1:2])$col covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss2_mix[,indi])} for (j in 11:40) #missing 1 parent pop1 {indi<-covindi(genokid[i,(j*2-1),1:2],genokid[i,(j*2),1:2])$col { if (sum(genokid[i,(j*2-1),3:4])==2) covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss1_mix_DD[,indi]) else if (sum(genokid[i,(j*2-1),3:4])==1) covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss1_mix_Dd[,indi]) else covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss1_mix_dd[,indi]) } } for (j in 41:50) #missing none {a<-covindi(genokid[i,(j*2-1),3:4],genokid[i,(j*2-1),5:6]) ##observed parental indi<-a$com row<-a$row indi0<-covindi(genokid[i,(j*2-1),1:2],genokid[i,(j*2),1:2])$col ##observed kids {if ((indi==1) || (indi==6)) ww<-as.matrix(z2kid[row,indi]) else ww<-diag(z2kid[row,indi])} temp<-as.matrix(x2kid[,row])%*%z2kid[row,indi]%*%solve(t(z2kid[row,indi])%*%solve(ww)%*%z2kid[row,indi])%*%t(z2kid[row,indi])%*%solve(ww) {if (indi==4) covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi0],temp[,indi]) else covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi0],temp[,1]) } } } p<-p2 q<-1-p w<-c(p^4,4*p^3*q,2*p^2*q^2,4*p^2*q^2,4*p*q^3,q^4) wmix<-diag(c(z2kid%*%w)) miss2_mix<-x2kid%*%z2kid%*%solve(t(z2kid)%*%solve(wmix)%*%z2kid)%*%t(z2kid)%*%solve(wmix) cDD<-c(1,2,4,5) xDD<-x2kid[,cDD]; z2DD<-z2kid[cDD,1:3] wDDv<-c(p^2,2*p*q,q^2); # Pr(g*|g**) wDD<-diag(c(z2DD%*%wDDv)) #sum_{g*}(Pr(g|g*)Pr(g*|g**)=Pr(g|g**) miss1_mix_DD0<-xDD%*%z2DD%*%solve(t(z2DD)%*%solve(wDD)%*%z2DD)%*%t(z2DD)%*%solve(wDD) miss1_mix_DD<-cbind(miss1_mix_DD0[,1:2],c(0,0),miss1_mix_DD0[,3:4],c(0,0),c(0,0),c(0,0),c(0,0)) wDd<-diag(c(z2kid[,c(2,4,5)]%*%wDDv)) miss1_mix_Dd<-x2kid%*%z2kid[,c(2,4,5)]%*%solve(t(z2kid[,c(2,4,5)])%*%solve(wDd)%*%z2kid[,c(2,4,5)])%*%t(z2kid[,c(2,4,5)])%*%solve(wDd) cdd<-c(5,6,8,9) xdd<-x2kid[,cdd]; z2dd<-z2kid[cdd,c(3,5,6)] wdd<-diag(c(z2dd%*%wDDv)) miss1_mix_dd0<-xdd%*%z2dd%*%solve(t(z2dd)%*%solve(wdd)%*%z2dd)%*%t(z2dd)%*%solve(wdd) miss1_mix_dd<-cbind(c(0,0),c(0,0),c(0,0),c(0,0),miss1_mix_dd0[,1:2],c(0,0),miss1_mix_dd0[,3:4]) for (i in 1:nrep) { for (j in 51:60) #missing 2 parents pop2 {indi<-covindi(genokid[i,(j*2-1),1:2],genokid[i,(j*2),1:2])$col covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss2_mix[,indi])} for (j in 61:90) #missing 1 parent pop2 {indi<-covindi(genokid[i,(j*2-1),1:2],genokid[i,(j*2),1:2])$col { if (sum(genokid[i,(j*2-1),3:4])==2) covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss1_mix_DD[,indi]) else if (sum(genokid[i,(j*2-1),3:4])==1) covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss1_mix_Dd[,indi]) else covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi],miss1_mix_dd[,indi]) } } for (j in 91:100) #missing none {a<-covindi(genokid[i,(j*2-1),3:4],genokid[i,(j*2-1),5:6]) ##observed parental indi<-a$com row<-a$row indi0<-covindi(genokid[i,(j*2-1),1:2],genokid[i,(j*2),1:2])$col ##observed kids {if ((indi==1) || (indi==6)) ww<-as.matrix(z2kid[row,indi]) else ww<-diag(z2kid[row,indi])} temp<-as.matrix(x2kid[,row])%*%z2kid[row,indi]%*%solve(t(z2kid[row,indi])%*%solve(ww)%*%z2kid[row,indi])%*%t(z2kid[row,indi])%*%solve(ww) {if (indi==4) covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi0],temp[,indi]) else covamis[i,(j*2-1):(j*2),]<-cbind(x2kid[,indi0],temp[,1]) } } } return(covamis) }