#This function takes data set. marker positions, omega, delta, #leta(log(1/lambda)), and locus (the location of x_s). #It will return (-(log-likelihood)) under alternative hypothesis ################################################################# #Other functions used inside of the function (they are attached at the end) #1.tr(cstage,nstage,dist,delta,lambda) #tr returns transition probability #from cstage to nstage #2.location(locus,pos) #Note: the length of pos is at least greater than or equal to two #If locus coincides with one of the marker positions, it will return #(locus,locus) #If locus is less than all other marker positions, #it will return (0,first marker) #If locus is greater than all other marker positions, #it will return (last marker,1) #If locus is between two markers, #it will return (left marker, right marker) ################################################################# #Basic algorithm #subtract all informative cases form data set #define: #rprob[i]=prob(L{i+1),L(i+1),,,L(m)|L(i)) #lprob[i]=prob(L(i-1),L(i-2),,,L(1)|L(i)) #compute rprob, lprob #If locus is one of the marker positions, #prob=omega^L(locus)*(1-omega)^(1-Locus)*rprob[locus]*lprob[locus] #Else #prob=omega*tr(1,L(right.marker))*rprob[right.marker] # *tr(1,L(left.marker))*lprob[left.marker] # +(1-omega)*tr(0,L(right.marker))*rprob[right.marker] # *tr(0,L(left.marker))*lprob[left.marker] #----------------------------------------------------- mloglik.alt<-function(delta,omega,leta,locus,data,pos){ lambda<-1/exp(leta) ntumor<-nrow(data) loglik<-0 for(i in 1:ntumor){ if((sum(data[i,]>=0))>=1) { #there exists at least one informative markers newpos<-pos[data[i,]>=0] #subtract informative markers newL<-data[i,data[i,]>=0] m<-length(newL) rprob<-rep(1,m) lprob<-rep(1,m) } else #there is no informative marker m<-0 #---- compute lprob and rprob------------------- if(m>=2) for(i in 1:(m-1)){ rprob[i]<-tr(newL[i],newL[i+1],(newpos[i+1]-newpos[i]),delta,lambda) lprob[i+1]<-tr(newL[i+1],newL[i],(newpos[i+1]-newpos[i]),delta,lambda) } if(m>1){ rprob<-rev(cumprod(rev(rprob))) lprob<-(cumprod(lprob)) # say, m=4. then # rprob<- c(p12*p23*p34*p45,p23*p34*p45,p34*p45, p45) # lprob <- c(p10,p10*p21,p10*p21*p32, p10*p21*p32*p43) } #--------------------------------------------------- if(m==1){ #only one informative marker if(locus==newpos[1]) #if locus coincides with newpos prob<-omega*newL[1]+(1-omega)*(1-newL[1]) #locus=marker position else { #if locus does not coincide with newpos dist<-abs(locus-newpos[1]) prob<-omega*tr(1,newL[1],dist,delta,lambda)+ (1-omega)*tr(0,newL[1],dist,delta,lambda) } } if(m>1){ #there are more than one informative markers temp<-location(locus,newpos) #get the location of two closest #informative markers(including 0,1) marker.left<-temp$left marker.right<-temp$right if(marker.left==marker.right){#locus=one of the informative markers index<-(newpos==locus) prob<-(omega*newL[index]*rprob[index]*lprob[index] +(1-omega)*(1-newL[index])*rprob[index]*lprob[index]) } else {#locus is different from all the marker location dist.right<-marker.right-locus dist.left<-locus-marker.left index.right<-(newpos==marker.right) index.left<-(newpos==marker.left) if(marker.left==0) #locus is the leftmost prob<-(omega*tr(1,newL[index.right],dist.right,delta,lambda) *rprob[index.right] +(1-omega)*tr(0,newL[index.right],dist.right,delta,lambda) *rprob[index.right]) else if(marker.right==1)#locus is the rightmost prob<-(omega*tr(1,newL[index.left],dist.left,delta,lambda) *lprob[index.left] +(1-omega)*tr(0,newL[index.left],dist.left,delta,lambda) *lprob[index.left]) else #locus is between two informative markers prob<-(omega*tr(1,newL[index.right],dist.right,delta,lambda) *rprob[index.right] *tr(1,newL[index.left],dist.left,delta,lambda)*lprob[index.left] +(1-omega)*tr(0,newL[index.right],dist.right,delta,lambda) *rprob[index.right] *tr(0,newL[index.left],dist.left,delta,lambda)*lprob[index.left]) } } if(m>=1) #to exclude the case where m=0 (no informative markers) loglik<-loglik+log(prob) } return(-loglik) #return minus loglik } ############################################################################# #This function takes L(x_1), L(x_2), distance between x1 and x2, delta, and #lambda #It returns transition probability tr<-function(cstage,nstage,dist,delta,lambda){ if(cstage==1){ if(nstage==1) return(1-(1-delta)*(1-exp(-lambda*dist))) else return((1-delta)*(1-exp(-lambda*dist))) } else { if(nstage==0) return(1-delta*(1-exp(-lambda*dist))) else return(delta*(1-exp(-lambda*dist))) } } ############################################################################# #This function returns two closet marker positions from locus #Note: length of pos is at least greater than or equal to two #If locus coincides with one of the marker positions, it will return #(locus,locus) #If locus is less than all other marker positions, #it will return (0,first marker) #If locus is greater than all other marker positions, #it will return (last marker,1) #If locus is between two markers, #it will return (left marker, right marker) location<-function(locus,pos){ if (sum(locus==pos)>=1) #if locus coincides with one of the marker positions return(list(left=locus,right=locus)) else{ aug.pos<-sort(c(0,pos,1,locus)) #add locus to marker position and sort index<-(locus==aug.pos) index.left<-c(index[-1],F) #move index by one to the left index.right<-c(F,index[-length(aug.pos)]) #move index by one to the right return(list(left=aug.pos[index.left],right=aug.pos[index.right])) } } #############################################################################