####################### function definition; adapted by Jason R Finley 2-8-2010 from Woods 2009 ######################## CCseJRFc <- function(X, Y, A, B, outputname) { #open as source file in R, then can use function. #on Mac OSX saves output to home directory of current user. #X and Y are two vectors (columns) of data for which we will compute a correlation #A and B are two other vectors for which we will compute a correlation. (note: one of A or B be could actually be the same vector as X or Y) #this function will also compute the standard error of the difference between the two tau-b correlations, both for the case in which XY and AB are independent and for cases in which they are dependent. the user will have to choose the appropriate values form the output. #there MUST be the same number of subjects in all 4 columns submitted to this function. #be sure to remove rows that have any empty cells before calling this function. n = length(X) ####### Declare variables for X and Y tijxy = matrix(c(0),nrow=n,ncol=n) #this is tau-sub-ijxy t2ijxy = matrix(c(0),nrow=n,ncol=n) #tau-sub-ijxy SQUARED tijxx = matrix(c(0),nrow=n,ncol=n) #this is tau-sub-ijxx, no ties on X tijyy = matrix(c(0),nrow=n,ncol=n) #this is tau-sub-ijyy, no ties on Y sumtidotxy = matrix(c(0),nrow=1,ncol=n) sumt2idotxy = matrix(c(0),nrow=1,ncol=n) sumtidotxx = matrix(c(0),nrow=1,ncol=n) sumtidotyy = matrix(c(0),nrow=1,ncol=n) tidotxy = matrix(c(0),nrow=1,ncol=n) t2idotxy = matrix(c(0),nrow=1,ncol=n) tidotxx = matrix(c(0),nrow=1,ncol=n) tidotyy = matrix(c(0),nrow=1,ncol=n) sumtijxy=sumt2ijxy=sumtijxx=sumtijyy=txx=tyy=txy=t2xy=0 vartidotxy=vart2idotxy=vartijxy=vart2ijxy=vartijxx=vartijyy=0 covXxXy=vartidotxx=varTxx=covYyXy=vartidotyy=varTyy=0 stixy=sti2xy=stijxy=st2ijxy=stixx=stijxx=stiyy=stijyy=0 varTxyC=varT2xyC=varTxxC=varTyyC=0 varbxyC=varGamxy=0 sebxyC=seGamxy=0 part1cxy=part2cxy=0 derivTxyGam=derivT2xyGam=0 txya=txyb=gammaxy=0 covijXxXy=covijYyXy=covij2XyXy=covidotXxXy=covidotYyXy=covidot2XyXy=covijXxXyC=covijYyXyC=covij2XyXyC=0 covidotXxXyC=covidotYyXyC=covidot2XyXyC=covXxXyC=covYyXyC=covXxYyC=cov2XyXyC=covidotXxYy=0 covijXxYy=covidotXxYyC=covijXxYyC=cov2XxXyC=covidot2XxXyC=covij2XxXyC=covidot2XxXy=0 covij2XxXy=cov2YyXyC=covidot2YyXyC=covij2YyXyC=covidot2YyXy=covij2YyXy=0 ##################### ####### Declare variables for A and B tijab = matrix(c(0),nrow=n,ncol=n) #this is tau-sub-ijab t2ijab = matrix(c(0),nrow=n,ncol=n) #tau-sub-ijab SQUARED tijaa = matrix(c(0),nrow=n,ncol=n) #this is tau-sub-ijaa, no ties on X tijbb = matrix(c(0),nrow=n,ncol=n) #this is tau-sub-ijbb, no ties on Y sumtidotab = matrix(c(0),nrow=1,ncol=n) sumt2idotab = matrix(c(0),nrow=1,ncol=n) sumtidotaa = matrix(c(0),nrow=1,ncol=n) sumtidotbb = matrix(c(0),nrow=1,ncol=n) tidotab = matrix(c(0),nrow=1,ncol=n) t2idotab = matrix(c(0),nrow=1,ncol=n) tidotaa = matrix(c(0),nrow=1,ncol=n) tidotbb = matrix(c(0),nrow=1,ncol=n) sumtijab=sumt2ijab=sumtijaa=sumtijbb=taa=tbb=tab=t2ab=0 vartidotab=vart2idotab=vartijab=vart2ijab=vartijaa=vartijbb=0 covAaAb=vartidotaa=varTaa=covBbAb=vartidotbb=varTbb=0 stiab=sti2ab=stijab=st2ijab=stiaa=stijaa=stibb=stijbb=0 varTabC=varT2abC=varTaaC=varTbbC=0 varbabC=varGamab=0 sebabC=seGamab=0 part1cab=part2cab=0 derivTabGam=derivT2abGam=0 taba=tabb=gammaab=0 covijAaAb=covijBbAb=covij2AbAb=covidotAaAb=covidotBbAb=covidot2AbAb=covijAaAbC=covijBbAbC=covij2AbAbC=0 covidotAaAbC=covidotBbAbC=covidot2AbAbC=covAaAbC=covBbAbC=covAaBbC=cov2AbAbC=covidotAaBb=0 covijAaBb=covidotAaBbC=covijAaBbC=cov2AaAbC=covidot2AaAbC=covij2AaAbC=covidot2AaAb=0 covij2AaAb=cov2BbAbC=covidot2BbAbC=covij2BbAbC=covidot2BbAb=covij2BbAb=0 ##################### #declare variables for covariances of X/Y and A/B covijXyAb = covijXyAa = covijXyBb = covijXxAb = covijXxAa = covijXxBb = covijYyAb = covijYyAa = covijYyBb = 0 covijXyAbC = covijXyAaC = covijXyBbC = covijXxAbC = covijXxAaC = covijXxBbC = covijYyAbC = covijYyAaC = covijYyBbC = 0 covidotXyAb = covidotXyAa = covidotXyBb = covidotXxAb = covidotXxAa = covidotXxBb = covidotYyAb = covidotYyAa = covidotYyBb = 0 covidotXyAbC = covidotXyAaC = covidotXyBbC = covidotXxAbC = covidotXxAaC = covidotXxBbC = covidotYyAbC = covidotYyAaC = covidotYyBbC = 0 covXyAbC = covXyAaC = covXyBbC = covXxAbC = covXxAaC = covXxBbC = covYyAbC = covYyAaC = covYyBbC = 0 ### Loop through all possible pairings for X & Y and count up concordances, discordances, ties for (i in 1:n) { for (j in 1:n) { #each ij pairing is only going to satisfy one of the following 5 IF statements... (I think) if ( ((X[i] > X[j]) && (Y[i] > Y[j])) || (X[i] < X[j]) && (Y[i] < Y[j])) { tijxy[i,j] = 1 #aka tij t2ijxy[i,j] = 1 #square of tij tijxx[i,j] = 1 tijyy[i,j] = 1 } #end if concordant if ( ((X[i] > X[j]) && (Y[i] < Y[j])) || (X[i] < X[j]) && (Y[i] > Y[j])) { tijxy[i,j] = -1 t2ijxy[i,j] = 1 tijxx[i,j] = 1 tijyy[i,j] = 1 } #end if discordant if ((X[i] == X[j]) && (Y[i] != Y[j])) { tijxy[i,j] = 0 t2ijxy[i,j] = 0 tijxx[i,j] = 0 tijyy[i,j] = 1 } #end if ties on X only if ((X[i] != X[j]) && (Y[i] == Y[j])) { tijxy[i,j] = 0 t2ijxy[i,j] = 0 tijxx[i,j] = 1 tijyy[i,j] = 0 } #end if ties on Y only if ((X[i] == X[j]) && (Y[i] == Y[j])) { tijxy[i,j] = 0 t2ijxy[i,j] = 0 tijxx[i,j] = 0 tijyy[i,j] = 0 } #end if ties on both X and Y } #end j } #end i ### Loop through all possible pairings for A & B and count up concordances, discordances, ties for (i in 1:n) { for (j in 1:n) { #each ij pairing is only going to satisfy one of the following 5 IF statements... (I think) if ( ((A[i] > A[j]) && (B[i] > B[j])) || (A[i] < A[j]) && (B[i] < B[j])) { tijab[i,j] = 1 #aka tij t2ijab[i,j] = 1 #square of tij tijaa[i,j] = 1 tijbb[i,j] = 1 } #end if concordant if ( ((A[i] > A[j]) && (B[i] < B[j])) || (A[i] < A[j]) && (B[i] > B[j])) { tijab[i,j] = -1 t2ijab[i,j] = 1 tijaa[i,j] = 1 tijbb[i,j] = 1 } #end if discordant if ((A[i] == A[j]) && (B[i] != B[j])) { tijab[i,j] = 0 t2ijab[i,j] = 0 tijaa[i,j] = 0 tijbb[i,j] = 1 } #end if ties on A only if ((A[i] != A[j]) && (B[i] == B[j])) { tijab[i,j] = 0 t2ijab[i,j] = 0 tijaa[i,j] = 1 tijbb[i,j] = 0 } #end if ties on B only if ((A[i] == A[j]) && (B[i] == B[j])) { tijab[i,j] = 0 t2ijab[i,j] = 0 tijaa[i,j] = 0 tijbb[i,j] = 0 } #end if ties on both A and B } #end j } #end i for (i in 1:n) { for (j in 1:n) { sumtijxy = sumtijxy + tijxy[i,j] sumt2ijxy = sumt2ijxy + t2ijxy[i,j] sumtijxx = sumtijxx + tijxx[i,j] sumtijyy = sumtijyy + tijyy[i,j] sumtidotxy[i] = sumtidotxy[i] + tijxy[i,j] sumt2idotxy[i] = sumt2idotxy[i] + t2ijxy[i,j] sumtidotxx[i] = sumtidotxx[i] + tijxx[i,j] sumtidotyy[i] = sumtidotyy[i] + tijyy[i,j] sumtijab = sumtijab + tijab[i,j] sumt2ijab = sumt2ijab + t2ijab[i,j] sumtijaa = sumtijaa + tijaa[i,j] sumtijbb = sumtijbb + tijbb[i,j] sumtidotab[i] = sumtidotab[i] + tijab[i,j] sumt2idotab[i] = sumt2idotab[i] + t2ijab[i,j] sumtidotaa[i] = sumtidotaa[i] + tijaa[i,j] sumtidotbb[i] = sumtidotbb[i] + tijbb[i,j] } } tidotxy = sumtidotxy/(n-1) #tidotxy for tau-xy = tau-a t2idotxy = sumt2idotxy/(n-1) #for tijxy squared tidotxx = sumtidotxx/(n-1) #tau sub i.xx for Somers' dyx tidotyy = sumtidotyy/(n-1) #for Somers' dxy txy = sumtijxy/(n*(n-1)) #this is tau-a t2xy = sumt2ijxy/(n*(n-1)) #(square tijxy then sum over pairs) txx = sumtijxx/(n*(n-1)) tyy = sumtijyy/(n*(n-1)) tidotab = sumtidotab/(n-1) #tidotab for tau-ab = tau-a t2idotab = sumt2idotab/(n-1) #for tijab squared tidotaa = sumtidotaa/(n-1) #tau sub i.aa for Somers' dba tidotbb = sumtidotbb/(n-1) #for Somers' dab tab = sumtijab/(n*(n-1)) #this is tau-a t2ab = sumt2ijab/(n*(n-1)) #(square tijab then sum over pairs) taa = sumtijaa/(n*(n-1)) tbb = sumtijbb/(n*(n-1)) #HERE1 covij #(co)variances of tau-ijxy, -ijxx, -ijyy for (i in 1:n) { for (j in 1:n) { if (i!=j) { #important for variances...covariances too stijxy = stijxy + ((tijxy[i,j] - txy)*(tijxy[i,j] - txy)) st2ijxy = st2ijxy + ((t2ijxy[i,j] - t2xy) *(t2ijxy[i,j] - t2xy)) stijxx = stijxx + ((tijxx[i,j] - txx)*(tijxx[i,j] - txx)) #for Somers' dyx stijyy = stijyy + ((tijyy[i,j] - tyy)*(tijyy[i,j] - tyy)) #for Somers' dxy covijXxXy = covijXxXy + ((tijxx[i,j] - txx)*(tijxy[i,j] - txy)) covijYyXy = covijYyXy + ((tijyy[i,j] - tyy)*(tijxy[i,j] - txy)) covijXxYy = covijXxYy + ((tijyy[i,j] - tyy)*(tijxx[i,j] - txx)) #for tau-b stijab = stijab + ((tijab[i,j] - tab)*(tijab[i,j] - tab)) st2ijab = st2ijab + ((t2ijab[i,j] - t2ab) *(t2ijab[i,j] - t2ab)) stijaa = stijaa + ((tijaa[i,j] - taa)*(tijaa[i,j] - taa)) #for Somers' dyx stijbb = stijbb + ((tijbb[i,j] - tbb)*(tijbb[i,j] - tbb)) #for Somers' dab covijAaAb = covijAaAb + ((tijaa[i,j] - taa)*(tijab[i,j] - tab)) covijBbAb = covijBbAb + ((tijbb[i,j] - tbb)*(tijab[i,j] - tab)) covijAaBb = covijAaBb + ((tijbb[i,j] - tbb)*(tijaa[i,j] - taa)) #for tau-b #components needed for covariances between X/Y and A/B covijXyAb = covijXyAb + ((tijxy[i,j] - txy)*(tijab[i,j] - tab)) covijXyAa = covijXyAa + ((tijxy[i,j] - txy)*(tijaa[i,j] - taa)) covijXyBb = covijXyBb + ((tijxy[i,j] - txy)*(tijbb[i,j] - tbb)) covijXxAb = covijXxAb + ((tijxx[i,j] - txx)*(tijab[i,j] - tab)) covijXxAa = covijXxAa + ((tijxx[i,j] - txx)*(tijaa[i,j] - taa)) covijXxBb = covijXxBb + ((tijxx[i,j] - txx)*(tijbb[i,j] - tbb)) covijYyAb = covijYyAb + ((tijyy[i,j] - tyy)*(tijab[i,j] - tab)) covijYyAa = covijYyAa + ((tijyy[i,j] - tyy)*(tijaa[i,j] - taa)) covijYyBb = covijYyBb + ((tijyy[i,j] - tyy)*(tijbb[i,j] - tbb)) } } } vartijxy = stijxy/((n*(n-1))-1) vart2ijxy = st2ijxy/((n*(n-1))-1) vartijxx = stijxx/((n*(n-1))-1) vartijyy = stijyy/((n*(n-1))-1) covijXxXyC = covijXxXy/((n*(n-1))-1) covijYyXyC = covijYyXy/((n*(n-1))-1) covijXxYyC = covijXxYy/((n*(n-1))-1) vartijab = stijab/((n*(n-1))-1) vart2ijab = st2ijab/((n*(n-1))-1) vartijaa = stijaa/((n*(n-1))-1) vartijbb = stijbb/((n*(n-1))-1) covijAaAbC = covijAaAb/((n*(n-1))-1) covijBbAbC = covijBbAb/((n*(n-1))-1) covijAaBbC = covijAaBb/((n*(n-1))-1) #components needed for covariances between X/Y and A/B covijXyAbC = covijXyAb/((n*(n-1))-1) covijXyAaC = covijXyAa/((n*(n-1))-1) covijXyBbC = covijXyBb/((n*(n-1))-1) covijXxAbC = covijXxAb/((n*(n-1))-1) covijXxAaC = covijXxAa/((n*(n-1))-1) covijXxBbC = covijXxBb/((n*(n-1))-1) covijYyAbC = covijYyAb/((n*(n-1))-1) covijYyAaC = covijYyAa/((n*(n-1))-1) covijYyBbC = covijYyBb/((n*(n-1))-1) #HERE2 covidot #(co)variances of tau-i.xy, -i.xx, -i.yy for (i in 1:n){ stixy = stixy + ((tidotxy[i] - txy) *(tidotxy[i] - txy)) sti2xy = sti2xy + ((t2idotxy[i] - t2xy) *(t2idotxy[i] - t2xy)) stixx = stixx + ((tidotxx[i] - txx)*(tidotxx[i] - txx)) #for Somers' dyx stiyy = stiyy + ((tidotyy[i] - tyy)*(tidotyy[i] - tyy)) #for Somers' dxy covidotXxXy = covidotXxXy + ((tidotxx[i] - txx)*(tidotxy[i] - txy)) covidotYyXy = covidotYyXy + ((tidotyy[i] - tyy)*(tidotxy[i] - txy)) covidotXxYy = covidotXxYy + ((tidotyy[i] - tyy)*(tidotxx[i] - txx)) #for tau-b stiab = stiab + ((tidotab[i] - tab) *(tidotab[i] - tab)) sti2ab = sti2ab + ((t2idotab[i] - t2ab) *(t2idotab[i] - t2ab)) stiaa = stiaa + ((tidotaa[i] - taa)*(tidotaa[i] - taa)) #for Somers' dba stibb = stibb + ((tidotbb[i] - tbb)*(tidotbb[i] - tbb)) #for Somers' dab covidotAaAb = covidotAaAb + ((tidotaa[i] - taa)*(tidotab[i] - tab)) covidotBbAb = covidotBbAb + ((tidotbb[i] - tbb)*(tidotab[i] - tab)) covidotAaBb = covidotAaBb + ((tidotbb[i] - tbb)*(tidotaa[i] - taa)) #for tau-b #components needed for covariances between X/Y and A/B covidotXyAb = covidotXyAb + ((tidotxy[i] - txy)*(tidotab[i] - tab)) covidotXyAa = covidotXyAa + ((tidotxy[i] - txy)*(tidotaa[i] - taa)) covidotXyBb = covidotXyBb + ((tidotxy[i] - txy)*(tidotbb[i] - tbb)) covidotXxAb = covidotXxAb + ((tidotxx[i] - txx)*(tidotab[i] - tab)) covidotXxAa = covidotXxAa + ((tidotxx[i] - txx)*(tidotaa[i] - taa)) covidotXxBb = covidotXxBb + ((tidotxx[i] - txx)*(tidotbb[i] - tbb)) covidotYyAb = covidotYyAb + ((tidotyy[i] - tyy)*(tidotab[i] - tab)) covidotYyAa = covidotYyAa + ((tidotyy[i] - tyy)*(tidotaa[i] - taa)) covidotYyBb = covidotYyBb + ((tidotyy[i] - tyy)*(tidotbb[i] - tbb)) } vartidotxy = stixy/(n-1) vart2idotxy = sti2xy/(n-1) #for tijxy squared vartidotxx = stixx/(n-1) vartidotyy = stiyy/(n-1) covidotXxXyC = covidotXxXy/(n-1) covidotYyXyC = covidotYyXy/(n-1) covidotXxYyC = covidotXxYy/(n-1) vartidotab = stiab/(n-1) vart2idotab = sti2ab/(n-1) #for tijab squared vartidotaa = stiaa/(n-1) vartidotbb = stibb/(n-1) covidotAaAbC = covidotAaAb/(n-1) covidotBbAbC = covidotBbAb/(n-1) covidotAaBbC = covidotAaBb/(n-1) #components needed for covariances between X/Y and A/B covidotXyAbC = covidotXyAb/(n-1) covidotXyAaC = covidotXyAa/(n-1) covidotXyBbC = covidotXyBb/(n-1) covidotXxAbC = covidotXxAb/(n-1) covidotXxAaC = covidotXxAa/(n-1) covidotXxBbC = covidotXxBb/(n-1) covidotYyAbC = covidotYyAb/(n-1) covidotYyAaC = covidotYyAa/(n-1) covidotYyBbC = covidotYyBb/(n-1) #CONSISTENT variances #for tau-a = tau-xy varTxyC = (((4*(n-2))*vartidotxy) + (2*vartijxy))/(n*(n-1)) seTxyC = sqrt(varTxyC) #no division by n #for tau-xx varTxxC = (((4*(n-2))*vartidotxx) + (2*vartijxx))/(n*(n-1)) #for tau-yy varTyyC = (((4*(n-2))*vartidotyy) + (2*vartijyy))/(n*(n-1)) #for tau-2 which is tau-xy squared varT2xyC = (((4*(n-2))*vart2idotxy) + (2*vart2ijxy))/(n*(n-1)) #for tau-a = tau-ab varTabC = (((4*(n-2))*vartidotab) + (2*vartijab))/(n*(n-1)) seTabC = sqrt(varTabC) #no division by n #for tau-aa varTaaC = (((4*(n-2))*vartidotaa) + (2*vartijaa))/(n*(n-1)) #for tau-bb varTbbC = (((4*(n-2))*vartidotbb) + (2*vartijbb))/(n*(n-1)) #for tau-2 which is tau-ab squared varT2abC = (((4*(n-2))*vart2idotab) + (2*vart2ijab))/(n*(n-1)) #CONSISTENT covariances #between tau-xx and tau-xy covXxXyC = ((4*(n-2)*covidotXxXyC) + (2*covijXxXyC))/(n*(n-1)) #between tau-yy and tau-xy covYyXyC = ((4*(n-2)*covidotYyXyC) + (2*covijYyXyC))/(n*(n-1)) #between tau-xx and tau-yy covXxYyC = ((4*(n-2)*covidotXxYyC) + (2*covijXxYyC))/(n*(n-1)) #between tau-aa and tau-ab covAaAbC = ((4*(n-2)*covidotAaAbC) + (2*covijAaAbC))/(n*(n-1)) #between tau-bb and tau-ab covBbAbC = ((4*(n-2)*covidotBbAbC) + (2*covijBbAbC))/(n*(n-1)) #between tau-aa and tau-bb covAaBbC = ((4*(n-2)*covidotAaBbC) + (2*covijAaBbC))/(n*(n-1)) #9 new covariances covXyAbC = ((4*(n-2)*covidotXyAbC) + (2*covijXyAbC))/(n*(n-1)) covXyAaC = ((4*(n-2)*covidotXyAaC) + (2*covijXyAaC))/(n*(n-1)) covXyBbC = ((4*(n-2)*covidotXyBbC) + (2*covijXyBbC))/(n*(n-1)) covXxAbC = ((4*(n-2)*covidotXxAbC) + (2*covijXxAbC))/(n*(n-1)) covXxAaC = ((4*(n-2)*covidotXxAaC) + (2*covijXxAaC))/(n*(n-1)) covXxBbC = ((4*(n-2)*covidotXxBbC) + (2*covijXxBbC))/(n*(n-1)) covYyAbC = ((4*(n-2)*covidotYyAbC) + (2*covijYyAbC))/(n*(n-1)) covYyAaC = ((4*(n-2)*covidotYyAaC) + (2*covijYyAaC))/(n*(n-1)) covYyBbC = ((4*(n-2)*covidotYyBbC) + (2*covijYyBbC))/(n*(n-1)) ## CC variances ## #gamma X Y derivTxyGam = 1/t2xy derivT2xyGam = -txy/(t2xy*t2xy) varGamxy = (derivTxyGam*derivTxyGam*varTxyC) + (derivT2xyGam*derivT2xyGam*varT2xyC) + (2*derivTxyGam*derivT2xyGam*cov2XyXyC) seGamxy = sqrt(varGamxy) #tau-b X Y #(from Cliff 1996 p.79 but includes denom. from Cliff & Charlin p. 700) part1cxy = varTxyC - (txy*((covXxXyC/txx)+(covYyXyC/tyy))) part2cxy = (txy*txy)*((varTxxC/(4*txx*txx))+(varTyyC/(4*tyy*tyy))+(covXxYyC/(2*txx*tyy))) varbxyC = (part1cxy + part2cxy)/(txx*tyy) sebxyC = sqrt(varbxyC) #point estimates of gamma family indices X Y taxy = txy tbxy = txy/(sqrt(txx*tyy)) gammaxy = txy/t2xy #gamma A B derivTabGam = 1/t2ab derivT2abGam = -tab/(t2ab*t2ab) varGamab = (derivTabGam*derivTabGam*varTabC) + (derivT2abGam*derivT2abGam*varT2abC) + (2*derivTabGam*derivT2abGam*cov2AbAbC) seGamab = sqrt(varGamab) #tau-b A B #(from Cliff 1996 p.79 but includes denom. from Cliff & Charlin p. 700) part1cab = varTabC - (tab*((covAaAbC/taa)+(covBbAbC/tbb))) part2cab = (tab*tab)*((varTaaC/(4*taa*taa))+(varTbbC/(4*tbb*tbb))+(covAaBbC/(2*taa*tbb))) varbabC = (part1cab + part2cab)/(taa*tbb) sebabC = sqrt(varbabC) #point estimates of gamma family indices A B taab = tab tbab = tab/(sqrt(taa*tbb)) gammaab = tab/t2ab #NOWWWW!!!! Get the covariance between the two tau-b coefficients!!!! #from Cliff & Charlin 1991, formula 20 mat1.1=(txx*tyy)^(-.5) mat1.2=(-.5)*(txy)*(txx^(-(3/2)))*(tyy^(-.5)) mat1.3=(-.5)*(txy)*(txx^(-.5))*(tyy^(-(3/2))) mat1 <- matrix(c(mat1.1,mat1.2,mat1.3),1,3) mat2 <- matrix(c(covXyAbC, covXyAaC, covXyBbC, covXxAbC, covXxAaC, covXxBbC, covYyAbC, covYyAaC, covYyBbC),3,3) mat3.1=(taa*tbb)^(-.5) mat3.2=(-.5)*(tab)*(taa^(-(3/2)))*(tbb^(-.5)) mat3.3=(-.5)*(tab)*(taa^(-.5))*(tbb^(-(3/2))) mat3 <- matrix(c(mat3.1,mat3.2,mat3.3),3,1) covarTaubXYAB <- mat1 %*% mat2 %*% mat3 seTaubXYABind <- sqrt((sebxyC^2) + (sebabC^2)) seTaubXYABdep <- sqrt(varbxyC + varbabC - (2*covarTaubXYAB)) zInd <- (tbxy-tbab)/seTaubXYABind zDep <- (tbxy-tbab)/seTaubXYABdep #Cliff 1996b, p 341, says to use z-test to compare taus. (other main option would be t-test, but I'm not sure what appropriate df would be. maybe: for ind: N1+N2-4. for dep with two taus drawing on same variable: N-3. for two dep with two taus drawing on different variables: N-4. pInd<-2*(1-pnorm(abs(zInd))) pDep<-2*(1-pnorm(abs(zDep))) outfile = file(description=paste(outputname,"_CompareTaub.txt",sep=""), open="a") #on Mac OSX saves to home directory of current user. #to output file cat( " gamma XY: ", gammaxy, " SE: ", seGamxy, "\n", "tau-a XY: ", taxy, " SE: ", seTxyC, "\n", "tau-b XY: ", tbxy, " SE: ", sebxyC, "\n", "gamma AB: ", gammaab, " SE: ", seGamab, "\n", "tau-a AB: ", taab, " SE: ", seTabC, "\n", "tau-b AB: ", tbab, " SE: ", sebabC, "\n\n", "SE (tau-b(XY) - tau-b(AB)) if independent: ", seTaubXYABind, "\n", "SE (tau-b(XY) - tau-b(AB)) if dependent: ", seTaubXYABdep, "\n", "indepdenent Z: ",zInd, " p: ",pInd,"\n", "dependent Z: ",zDep, " p: ",pDep,"\n", file=outfile,append=TRUE) close(outfile) } #end