1
$\begingroup$

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:

enter image description here

enter image description here

(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?

Thanks in advance.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Glen_b -Reinstate Monica Oct 8 '13 at 4:38
  • $\begingroup$ Using ASE for this data is not my idea :) I am using the data presented in the article mentioned in my question. $\endgroup$ – mjaniec Oct 8 '13 at 4:44
  • $\begingroup$ @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. $\endgroup$ – mjaniec Oct 8 '13 at 5:10
  • $\begingroup$ 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. $\endgroup$ – ttnphns Oct 8 '13 at 7:57
  • $\begingroup$ 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. $\endgroup$ – Frank Harrell Oct 8 '13 at 11:43
0
$\begingroup$

After some additional research, I was able to implement Somers' D ASE in R and, additionally, in F#.

Both codes are not the most clean, but they seem to work fine:

http://www.reakkt.com/2013/12/sommersd-two-dirty-implementations-r-vs.html

I would like to thank all who assisted me in my search for the solution :)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.