version 16.0 mata: real colvector tidottree(real colvector x, real colvector y, | real colvector weight, real colvector xcen, real colvector ycen) { /* Return weighted concordance-discordance difference counts between an X-variable and a Y-variable (both possibly censored), using a search tree algorithm (WHICH ASSUMES THE DATA ARE SORTED BY X). x is the X-variable. y is the Y-variable. weight contains the weights used to calculate tidot. xcen contains the censorship status of the X-variable (-1 for left-censored, 0 for uncensored, 1 for right-censored). ycen contains the censorship status of the Y-variable (-1 for left-censored, 0 for uncensored, 1 for right-censored). *! Author: Roger Newson *! Date: 15 September 2018 */ real scalar narg, nobs, nxval, nyval, i1, i2, i2min, i2max, root, ycur, xcencur, ycencur, wtcur, llcur, lgcur, glcur, ggcur, nodecur, yvalcur real matrix xpanel, ytree real colvector tidot, yval, swlt, sweqlc, sweqnc, sweqrc, swgt /* Set number of observations */ nobs=rows(x) /* Initialise absent arguments. */ narg=args(); if (narg<3) {;weight=J(nobs,1,1);} if (narg<4) {;xcen=J(nobs,1,0);} if (narg<5) {;ycen=J(nobs,1,0);} /* Return empty output in event of empty input. */ if(nobs<1) { return(J(0,1,0)) } /* Initialize X-panel setup matrix and Y-value tree. */ xpanel=panelsetup(x,1) nxval=rows(xpanel) yval=uniqrows(y) nyval=rows(yval) ytree=blncdtree(nyval) root=trunc((1+nyval)/2) /* Initialize tidot and node sums of weights to zero. */ tidot=J(nobs,1,0) swlt=J(nyval,1,0) sweqlc=J(nyval,1,0) sweqnc=J(nyval,1,0) sweqrc=J(nyval,1,0) swgt=J(nyval,1,0) /* Iterate over X-values from the lowest to the highest, extracting lower subtotals llcur and lgcur from tree to update tidot and then updating lower tree node sums of weights using weight. */ for(i1=1;i1<=nxval;i1++) { /* Beginning of iteration for current X-value. */ i2min=xpanel[i1,1] i2max=xpanel[i1,2] /* First observation loop for current X-value. Extract lower subtotals from tree and update tidot. */ for(i2=i2min;i2<=i2max;i2++) { /* Beginning of iteration for current observation. */ ycur=y[i2] xcencur=xcen[i2] ycencur=ycen[i2] if (xcencur>=0) { /* Extract left and right subtotals from tree and update tidot for current observation. */ llcur=0 lgcur=0 nodecur=root while (nodecur>0) { yvalcur=yval[nodecur] if (ycuryvalcur) { /* Enter right subtree */ if (ycencur>=0) llcur = quadsum((llcur,swlt[nodecur],(sweqlc[nodecur]+sweqnc[nodecur]))) nodecur=ytree[nodecur,2] } else { /* Current node has current Y-value */ if (ycencur<=0) lgcur = quadsum((lgcur,swgt[nodecur])) if (ycencur>=0) llcur = quadsum((llcur,swlt[nodecur])) nodecur=0 } } tidot[i2] = quadsum((tidot[i2],(llcur-lgcur))) } /* End of iteration for current observation. */ } /* Second observation loop for current X-value. Update lower node sums of weights in tree with weights in weight. */ for(i2=i2min;i2<=i2max;i2++) { /* Beginning of iteration for current observation. */ wtcur=weight[i2] ycur=y[i2] xcencur=xcen[i2] ycencur=ycen[i2] if ((!missing(wtcur)) & (xcencur<=0)) { /* Update lower node sums of weights from tree with current observation. */ nodecur=root while(nodecur>0) { yvalcur=yval[nodecur] if (ycuryvalcur) { /* Enter right subtree */ if (ycencur>=0) swgt[nodecur] = quadsum((swgt[nodecur],wtcur)) nodecur=ytree[nodecur,2] } else { /* Current node has current Y-value */ if(ycencur<0) sweqlc[nodecur] = quadsum((sweqlc[nodecur],wtcur)) else if (ycencur>0) sweqrc[nodecur] = quadsum((sweqrc[nodecur],wtcur)) else sweqnc[nodecur] = quadsum((sweqnc[nodecur],wtcur)) nodecur=0 } } /* End of iteration for current observation. */ } } /* End of iteration for current X-value. */ } /* Re-initialize node sums of weights to zero. */ swlt=J(nyval,1,0) sweqlc=J(nyval,1,0) sweqnc=J(nyval,1,0) sweqrc=J(nyval,1,0) swgt=J(nyval,1,0) /* Iterate over X-values from the highest to the lowest, extracting upper subtotals glcur and ggcur from tree to update tidot and then updating upper tree node sums of weights using weight. */ for(i1=nxval;i1>=1;i1--) { /* Beginning of iteration for current X-value. */ i2min=xpanel[i1,1] i2max=xpanel[i1,2] /* First observation loop for current X-value. Extract upper subtotals from tree and update tidot. */ for(i2=i2max;i2>=i2min;i2--) { /* Beginning of iteration for current observation. */ ycur=y[i2] xcencur=xcen[i2] ycencur=ycen[i2] if (xcencur<=0) { /* Extract left and right subtotals from tree and update tidot for current observation. */ glcur=0 ggcur=0 nodecur=root while (nodecur>0) { yvalcur=yval[nodecur] if (ycuryvalcur) { /* Enter right subtree */ if (ycencur>=0) glcur = quadsum((glcur,swlt[nodecur],(sweqlc[nodecur]+sweqnc[nodecur]))) nodecur=ytree[nodecur,2] } else { /* Current node has current Y-value */ if (ycencur<=0) ggcur = quadsum((ggcur,swgt[nodecur])) if (ycencur>=0) glcur = quadsum((glcur,swlt[nodecur])) nodecur=0 } } tidot[i2] = tidot[i2] + (ggcur-glcur) } /* End of iteration for current observation. */ } /* Second observation loop for current X-value. Update lower node sums of weights in tree with weights in weight. */ for(i2=i2max;i2>=i2min;i2--) { /* Beginning of iteration for current observation. */ wtcur=weight[i2] ycur=y[i2] xcencur=xcen[i2] ycencur=ycen[i2] if ((!missing(wtcur)) & (xcencur>=0)) { /* Update upper node sums of weights from tree with current observation. */ nodecur=root while(nodecur>0) { yvalcur=yval[nodecur] if (ycuryvalcur) { /* Enter right subtree */ if (ycencur>=0) swgt[nodecur] = quadsum((swgt[nodecur],wtcur)) nodecur=ytree[nodecur,2] } else { /* Current node has current Y-value */ if(ycencur<0) sweqlc[nodecur] = quadsum((sweqlc[nodecur],wtcur)) else if (ycencur>0) sweqrc[nodecur] = quadsum((sweqrc[nodecur],wtcur)) else sweqnc[nodecur] = quadsum((sweqnc[nodecur],wtcur)) nodecur=0 } } } /* End of iteration for current observation. */ } /* End of iteration for current X-value. */ } /* Return tidot as result */ return(tidot) } end