R/wdbm.R

In pmr: Probability Models for Ranking Data

wdbm <- function(dset,dtype="tau"){

nitem <- ncol(dset)-1

test <- matrix(data = 0, nrow = factorial(nitem), ncol = nitem, byrow = TRUE)
temp1 <- 1:nitem
i <- 1
w <- rep(1,nitem)
modal <- 1:nitem

## generate a list of all possible rankings
for (j in 1:(nitem^nitem-1)){
temp1[1] <- nitem - j%%nitem
temp2 <- j - j%%nitem
for (k in nitem:2){
temp1[k] <- nitem - temp2%/%(nitem^(k-1))
temp2 <- temp2 - (nitem-temp1[k])*(nitem^(k-1))
}
temp2 <- 0
for (l in 1:nitem){
for (m in 1:nitem){
if (temp1[l] == temp1[m] && l != m){
temp2 <- 1
}
}
}
if (temp2 == 0){
for (p in 1:nitem){
test[i,p] = temp1[p]
}
i <- i+1
}
}

n <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
for (k in 1:nrow(dset)){
temp_ind <- 0
for (l in 1:nitem){
if (test[j,l] != dset[k,l]) {temp_ind <- temp_ind + 1}
}
if (temp_ind == 0) {n[j] <- dset[k,nitem+1]}
}
}
test2 <- cbind(test, n)

## define distance measure
if (dtype == "rho"){
wdist <- function(x,y,w){
d <- 0
for (j in 1:nitem){
d <- d + (x[j] - y[j])^2 *w[modal[j]]
}
d^0.5
}
}

if (dtype == "rho2"){
wdist <- function(x,y,w){
d <- 0
for (j in 1:nitem){
d <- d + (x[j] - y[j])^2 *w[modal[j]]
}
d
}
}

if (dtype == "foot"){
wdist <- function(x,y,w){
d <- 0
for (j in 1:nitem){
d <- d + abs(x[j] - y[j]) *as.numeric(w[modal[j]])
}
d
}
}

if (dtype == "tau"){
wdist <- function(x,y,w){
d <- 0
for (j in 1:(nitem-1)){
for (k in (j+1):nitem){
if ((x[j]-x[k])*(y[j]-y[k]) < 0) {d <- d+w[modal[k]]*w[modal[j]]}
}
}
d
}
}

tempw <- rep(1,nitem)
if (dtype=="phicom"){tempw <- rep(1,(nitem-1))}
td <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
td[j] <- 0
for (k in 1:factorial(nitem)){
td[j] <- td[j] + n[k]*wdist(test[j,],test[k,],tempw)
}
}
test3 <- cbind(test2, td)

## compute temporary modal ranking
temp1 <- max(td)
for (j in 1:factorial(nitem)){
if (td[j] == min(td)){
for (k in 1:nitem){
modal[k] <- test3[j,k]
}
}
}

temp_modal <- 1:nitem
mup <- 1
while (mup == 1){
## loglikelihood function

loglik_wdbm <- function(lambda){
ed <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
ed[j] <- exp(-wdist(modal,test3[j,1:nitem],lambda))
}
pc <- sum(ed)
pr <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
pr[j] <- ed[j]/pc
}
ll <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
ll[j] <- -log(pr[j])*test3[j,(nitem+1)]
}
sum(ll)
}
if (dtype!="phicom") {up <- optim(rep(1,nitem), loglik_wdbm, NULL, method = "BFGS", hessian = TRUE)}
## if (dtype!="phicom") {up <- mle(minuslogl = loglik_wdbm, start = list(lambda=rep(1,nitem)), method = "BFGS", nobs=as.integer(sum(dset[,nitem+1])))}
w_up <- up$par
for (j in 1:nitem){
temp_modal[j] <- modal[j]
}

twd <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
twd[j] <- 0
for (k in 1:factorial(nitem)){
twd[j] <- twd[j] + n[k]*wdist(test[j,],test[k,],w_up)
}
}
test3 <- cbind(test2, twd)

temp1 <- max(td)
for (j in 1:factorial(nitem)){
if (twd[j] == min(twd)){
for (k in 1:nitem){
modal[k] <- test3[j,k]
}
}
}
mup <- 0
for (j in 1:nitem){
if (temp_modal[j] != modal[j]) {mup <- 1}
}
}

## compute expected value
ed <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
ed[j] <- exp(-wdist(modal,test3[j,1:nitem],up$par))
}
pc <- sum(ed)
fitted <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
fitted[j] <- ed[j]/pc * sum(test3[,(nitem+1)])
}
ss <- rep(0,factorial(nitem))
for (j in 1:factorial(nitem)){
ss[j] <- (fitted[j] - test3[j,(nitem+1)])^2/fitted[j]
}

#lst <- list(modal.ranking=modal, loglik=up$value, par=up$par, se=(diag(solve(up$hessian)))^0.5, fit.value=fitted, residual=sum(ss))
#return(lst)
##require(stats4)
message("Maximum Likelihood Estimation of the Weighted Distance-based Model")
dtype_full <- "Weighted Kendall's tau"
if (dtype == "rho") dtype_full <- "Weighted Spearman's rho"
if (dtype == "rho2") dtype_full <- "Weighted Spearman's rho square"
if (dtype == "foot") dtype_full <- "Weighted Spearman's footrule"
message("Distance type: ", dtype_full)
message("Modal ranking: ", modal)
message("Chi-square residual statistic: ", round(sum(ss), digits = 2), ", df: ", factorial(nitem))
out2 <- new("mle")
out2@coef <- up$par
out2@fullcoef <- up$par
out2@vcov <- solve(up$hessian)
out2@min <- up$value
out2@details <- up
out2@minuslogl <- loglik_wdbm
out2@method <- "BFGS"
return(out2)
}