# Problem with calculating Asymptotic Standard Error for Somer's D

I'm trying to calculate Asymptotic Standard Error (ASE) for Somer's D measure of association for some population data.

The formula for variation (sd^2) I've found (in the article specified below) is:

(this is identical with equations presented in "Analysis of Ordinal Categorical Data" by Alan Agresti)

I've tried to implement it in R:

pairsC <- function(x,prob=FALSE) {

l <- nrow(x) # rows
k <- ncol(x) # columns

pc <- 0

for (i in 1:l)
for (j in 1:k) {

row.idx <- (1:l)[which((1:l)>i)]
col.idx <- (1:k)[which((1:k)>j)]

if (length(row.idx)>0 && length(col.idx)>0)

pc <- pc+x[i,j]*sum(x[row.idx,col.idx])

}

ifelse(prob,2*pc,pc)

}

pairsD <- function(x,prob=FALSE) {

l <- nrow(x) # rows
k <- ncol(x) # columns

pd <- 0

for (i in 1:l)
for (j in 1:k) {

row.idx <- (1:l)[which((1:l)>i)]
col.idx <- (1:k)[which((1:k)<j)]

if (length(row.idx)>0 && length(col.idx)>0)

pd <- pd+x[i,j]*sum(x[row.idx,col.idx])

}

ifelse(prob,2*pd,pd)

}

prob <- TRUE

pc <- function(x,i,j) {

pcx <- 0

l <- nrow(x) # rows
k <- ncol(x) # columns

for a in 1:l
for b in 1:k

if (a<i and b<j) pcx <- pcx+x[a,b]

for a in 1:l
for b in 1:k

if (a>i and b>j) pcx <- pcx+x[a,b]

pcx

}

pd <- function(x,i,j) {

pdx <- 0

l <- nrow(x) # rows
k <- ncol(x) # columns

for a in 1:l
for b in 1:k

if (a<i and b>j) pdx <- pdx+x[a,b]

for a in 1:l
for b in 1:k

if (a>i and b<j) pdx <- pdx+x[a,b]

pdx

}

fi <- function(x,i,j) {

l <- nrow(x) # rows
k <- ncol(x) # columns

tmp <-

-2 * sum(x[i,]) * (pairsC(x,prob)-pairsD(x,prob)) +

-2 * (1-sum(x[setdiff(1:l,i),])^2) * (pc(x,i,j)-pd(x,i,j))

ifelse(length(tmp>0),tmp,0)

}

# nom1 & nom2 - nominators

l <- nrow(x) # rows
k <- ncol(x) # columns

nom1 <- 0
nom2 <- 0

for (i in 1:l)
for (j in 1:k) {

nom1 <- nom1 + x[i,j]*fi(x,i,j)^2

nom2 <- nom2 + x[i,j]*fi(x,i,j)

}

nom2 <- nom2^2

# den - denominator

den <- (1-sum(rowSums(x)^2))^2

# ASE

sqrt(nom1-nom2)/den


The table I'm trying to calculate ASE for is:

> portfolio
[,1]   [,2]  [,3]   [,4]  [,5]   [,6]   [,7]
defaulted    0.000025 0.0001 0.001 0.0025 0.004 0.0075 0.0125
nondefaulted 0.049975 0.0999 0.199 0.2475 0.196 0.1425 0.0375


According to an article I've taken the table from (Oliver Blumke "A proposal for a validation methodology for the discriminatory power of a rating system over time", The Journal of Risk Model Validation, Vol. 5/Num. 1, Spring 2011), the Somer's D should be 0.684 and ASE 2.378.

While Somer's D I've calculated is identical to the one in the article, less the sign, i.e. -0.684 instead of +0.684, "my" ASE is 7.5468.

Can anyone help me with deriving the correct ASE value?

• The difference in sign may just be an interchange of x and y; I don't think that in any way explains the difference in ASE though. ... why would you use ASE for a sample size of 7 though? – Glen_b -Reinstate Monica Oct 8 '13 at 4:38
• Using ASE for this data is not my idea :) I am using the data presented in the article mentioned in my question. – mjaniec Oct 8 '13 at 4:44
• @Glen_b You are welcome :) BTW: you are right about the sign of the measure and interchange of some data. When I use -pairsC(x,prob)+pairsD(x,prob) instead of pairsC(x,prob)-pairsD(x,prob), I get +0.684. – mjaniec Oct 8 '13 at 5:10
• I've checked your data in SPSS. (Since you didn't indicate the sample size N, I took N=1, i.e. your proportion values are taken as frequences.) Somer's D (columns dependent) = -.684; ASE (not assuming H0) = 2.378. – ttnphns Oct 8 '13 at 7:57
• The R Hmisc package's rcorr.cens function uses the standard $U$-statistic standard error estimator for $D_{xy}$, calculated efficiently using Fortran. You might take a look. – Frank Harrell Oct 8 '13 at 11:43