From: Wuming Gong <wuming.gong_at_gmail.com>

Date: Fri 02 Sep 2005 - 16:49:47 EST

}

gamma <- (C - D) / (C + D) # gamma = 0.221 for this example.

}

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Sep 02 17:02:14 2005

Date: Fri 02 Sep 2005 - 16:49:47 EST

Dear list,

I have a problem on calculating the standard error of Goodman-Kurskal's gamma using delta method. I exactly follow the method and forumla described in Problem 3.27 of Alan Agresti's Categorical Data Analysis (2nd edition). The data I used is also from the job satisfaction vs. income example from that book.

job <- matrix(c(1, 3, 10, 6, 2, 3, 10, 7, 1, 6, 14, 12, 0, 1, 9, 11), nrow = 4, ncol = 4, byrow = TRUE, dimnames = list(c("< 15,000", "15,000 - 25,000", "25,000 - 40,000", "> 40,000"), c("VD", "LD", "MS","VS")))

The following code is for calculating gamma value, which is consistent with the result presented in section 2.4.5 of that book.

C <- 0

D <- 0

for (i in 1:nrow(job)){

for (j in 1:ncol(job)){ pi_c <- 0 pi_d <- 0 for (h in 1:nrow(job)){ for (k in 1:ncol(job)){ if ((h > i & k > j) | (h < i & k < j)){ pi_c <- pi_c + job[h, k]/sum(job)

}

if ((h > i & k < j) | (h < i & k > j)){ pi_d <- pi_d + job[h, k]/sum(job)

}

} } C <- C + job[i, j] * pi_c D <- D + job[i, j] * pi_d }

}

gamma <- (C - D) / (C + D) # gamma = 0.221 for this example.

The following code is for calculating stardard error of gamma.
sigma.squared <- 0

for (i in 1:nrow(job)){

for (j in 1:ncol(job)){ pi_c <- 0 pi_d <- 0 for (h in 1:nrow(job)){ for (k in 1:ncol(job)){ if ((h > i & k > j) | (h < i & k < j)){ pi_c <- pi_c + job[h, k]/sum(job)

}

if ((h > i & k < j) | (h < i & k > j)){ pi_d <- pi_d + job[h, k]/sum(job)

}

} } phi <- 4 * (pi_c * D - pi_d * C) / (C + D)^2 sigma.squared <- sigma.squared + phi^2 }

}

se <- (sigma.squared/sum(job))^.5 # 0.00748, which is different from the SE 0.117 given in section 3.4.3 of that book.

I am not able to figure out what is the problem with my code... Could anyone point out what the problem is?

Thanks.

Wuming

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Sep 02 17:02:14 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:40:02 EST
*