Upload New File

module4/Reproduction/OPRC1.R

+#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+# Author of the R codes: Jen-Chieh Teng.
+# For the paper submitted to CSDA: Lim, A., Chiang, C. T., Teng, J.C., 2020. Estimating robot strengths with application to selection of alliance members in FIRST robotics competitions. (Under Revision)
+# 
+#
+# These are NOT but a user-friendly version of those used in the paper. 
+# Please inform me of any errors that I might have introduced in trying to make the codes more readable. 
+# These codes are self-contained in such a way that parameter estimates are readily available once the 
+# user loads their dataset formatted in the appropriate structure.
+#---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+library(readxl)
+
+#################################################################################################
+# the functions of robot strength estimates and estimated rank correlation of two variables     #
+#################################################################################################
+
+bhat.fn = function( x , y ){
+  xc = x %*% INF1   
+  bhatc = solve( t(xc) %*% xc ) %*% t(xc) %*% y
+  binf = INF1 %*% bhatc
+  binf}
+
+RC = function( x1 , x2 ){
+  t = length(x1)
+  r1 = outer(x1,x1,FUN="-")
+  r2 = outer(x2,x2,FUN="-")
+  RC = (sum(sign(r1*r2)>0)  + 0.5*(sum(sign(r1*r2)==0) -t ))  / (t*(t-1))
+  RC
+}
+
+##################################################################################################
+#     Import excel files of a specific division of a tournament in the qualification stage (e.g. #        
+#     Archimedes in 2018_Detroit.xlsx or Carver in 2018_Houston.xlsx) and the playoff stage      #       
+#     (e.g. Archimedes in 2018_Detroit_Playoffs.xlsx or Carver in 2018_Houston_Playoffs.xlsx)    #                                      
+##################################################################################################
+
+  i = 1 # substitute i=1 by i=2,...,or i=6 for the rest divisions 
+  Qualification <- read_excel("2018-Detroit.xlsx", sheet = i) 
+  Playoffs <- read_excel("2018-Detroit-Playoffs.xlsx" , sheet = i)
+  # The input files are named 2018-Detroit.xlsx and Data_2018-Detroit-Playoffs.xlsx (or 2018-Houston.xlsx 
+  # and Data_2018-Houston-Playoffs.xlsx) with the following format:
+  # the 1st column: match numbers
+  # the 2nd to 4th columns: robot IDs in the red alliance
+  # the 5th to 7th columns: robot IDs in the blue alliance
+  # the 8th column: a score of the red alliace
+  # the 9th column: a score of the blue alliance
+  # the 10th column: an adjusted score of the red alliance
+  # the 11th column: an adjusted score of the blue allaince
+  # the 15th column: match numbers
+  # the remaining columns: covariate matrix X with X_ij being equal to 1 if the jth robot at the ith match 
+  # assigned to the red alliance; -1 if the jth robot at the ith match assigned to the blue alliance; and 
+  # 0 otherwise. 
+  
+  
+  # rename 2018_Detroit.xlsx (or 2018_Houston.xlsx) as Data 
+  Data = data.frame(Qualification)  
+  
+  # rename 2018_Detroit_Playoffs.xlsx (or 2018_Houston_Playoffs.xlsx) as Data 
+  Data_po = data.frame(Playoffs)  
+  
+  # number of matches (M), number of robotics teams (K), and IDs in the qualification stage 
+  M = length(Data[,1]) 
+  K = length(Data[1,]) - 15
+  ID0 = colnames(Qualification)[16:(K+15)]
+  ID = cbind(as.numeric(ID0),seq(1,K))
+  
+  # the minimum number (m0) of matches played by the robotics teams
+  m0 = ceiling(M*6/K)-1
+  
+  # the corresponding response y0 and covariates x0, which are derived from qualification data
+  yr = as.matrix(Data[1:M,8])
+  yb = as.matrix(Data[1:M,9])
+  y0 = rbind(yr,yb)
+  ny = length(y0)
+  x00 = as.matrix(Data[1:M,1:K+15])
+  xr = ( x00 == 1 )*1
+  xb = ( x00 == -1 )*1
+  x0 = rbind(xr,xb)
+  
+  # the number of matches (Mpo) and number of robotics teams (K0po) in the playoff stage
+  Mpo = length(t(Data_po[,1])) 
+  K0po = length(Data_po[1,]) - 15
+  
+  # the corresponding response y0po, covariates x0po, and IDs, which are derived from playoff data 
+  yrpo = as.matrix(Data_po[1:Mpo,8])
+  ybpo = as.matrix(Data_po[1:Mpo,9])
+  y0po = rbind(yrpo,ybpo)
+  x0po0 = as.matrix(Data_po[1:Mpo,1:K0po+15])
+  xrpo = ( x0po0 == 1 )*1
+  xbpo = ( x0po0 == -1 )*1
+  x0po = rbind(xrpo,xbpo)
+  nypo = length(y0po)
+  dim(y0po) = c(nypo,1)
+  Kpo = sum((colSums(abs(x0po)) > 0)*1)
+  po = colSums((abs(x0po)))
+  dim(po) = c(K,1)
+  poinf = cbind(po,ID)
+  row_sub = apply ( poinf, 1 ,function(row) all(row !=0 ))
+  poat = as.matrix(poinf[row_sub, ])
+  ID.po = poat[,2:3]
+  
+  ##################################################################################################
+  # the LSE of robot strengths (bhat0) in the OPR model                                            #
+  ##################################################################################################
+  bhat0 = solve( t(x0) %*% x0 ) %*% t(x0) %*% y0
+  
+  ##################################################################################################
+  # the information matrices for clusters (cluster.inf0) and parameter estimates (bhat.inf0)       #
+  # under variant cluster numbers = 2,…,K (clus_n) for the proposed methods 1                      #      
+  ################################################################################################## 
+  clus_n = seq(2,K)
+  l_n = length( clus_n )
+  cluster.inf0 = bhat.inf0 =  matrix( 0 , K , l_n)
+  fit = hclust( dist(bhat0) , method = "centroid" )
+  cluster.inf0[1:K,l_n] = cutree( fit , k = clus_n[l_n])
+  bhat.inf0[1:K,l_n] = bhat0
+  
+  for(w in 1:(l_n -1)){
+    grp = cutree( fit , k = clus_n[(l_n-w)])
+    grp0 = sort(unique(grp))
+    lgrp0 = length(grp0)
+    INF1 = outer( grp , grp0 , FUN="==") *1 
+    cluster.inf0[1:K , ( l_n-w)] = grp
+    bhat.inf0[ 1:K, ( l_n-w)] = bhat.fn(x0,y0)
+  }
+  
+  
+  
+  ##################################################################################################
+  # predictive matrices PR’s and MSPE’s, which are computed based on qualification data,           #
+  # under variant combinations of cluster number (c) and number of matches (M_l), l = 6,…,m0       #
+  ##################################################################################################
+  m = seq(6,m0)
+  l_m = length(m)
+  c0 = ceiling(m * K /6)
+  c0 = c(c0 , c0[l_m] +1 )
+  
+  ##################################################################################################
+  #                                     Method 1                                                   #                                   
+  ##################################################################################################
+    PRhatcv = MSPEcv = matrix (0 , l_n , l_m)
+    
+    for( w1 in 1:l_m){ #crossvalidation
+      yrm = yr[-c( (c0[w1] +1) : c0[l_m+1] )]
+      ybm = yb[-c( (c0[w1] +1) : c0[l_m+1] )]
+      nym = length(yrm) 
+      xrm = xr[-c(c0[w1] + 1 : c0[l_m+1]) ,]
+      xbm = xb[-c(c0[w1] + 1 : c0[l_m+1]) ,]
+      
+      for(w2 in 1:l_n){
+        grp = cluster.inf0[(1:K),w2]
+        grp0 = sort(unique(grp))
+        lgrp0 = length(grp0)
+        INF1 = outer( grp , grp0 , FUN="==") *1
+        
+        # estimated PR and MSPE
+        
+        ycv = c(ybm, yrm)
+        nycvm = length(ycv)
+        dim(ycv) = c(nycvm,1)
+        dy_cv = yrm-ybm
+        
+        xcv = rbind(xbm,xrm)%*%INF1
+        H0=solve(t(xcv)%*%xcv)%*%t(xcv)
+        H_cv=xcv%*%H0
+        
+        bhat_cv=H0%*%ycv
+        pr_cv0=(xcv[-c(1:nym),]-xcv[1:nym,])%*%bhat_cv
+        pr_cv0=pr_cv0[1:nym]
+        resid_cv0=ycv-H_cv%*%ycv
+        
+        w01=diag(H_cv)
+        w02=diag(H_cv[1:nym,-c(1:nym)])
+        wd_cv=(1-w01[1:nym])*(1-w01[-c(1:nym)])-w02^2
+        w1_cv=((1-w01[-c(1:nym)])*resid_cv0[1:nym]+w02*resid_cv0[-c(1:nym)])/wd_cv
+        w2_cv=((1-w01[1:nym])*resid_cv0[-c(1:nym)]+w02*resid_cv0[1:nym])/wd_cv
+        pr_cv=pr_cv0-t(t(H_cv[-c(1:nym),1:nym]-H_cv[1:nym,1:nym])*w1_cv)-t(t(H_cv[-c(1:nym),-c(1:nym)]-H_cv[1:nym,-c(1:nym)])*w2_cv)
+        
+        
+        
+        resid_cv=dy_cv-pr_cv
+        Fhat_cv0=(t(t(resid_cv)+diag(pr_cv))<=0)*1
+        diag(Fhat_cv0)=0
+        Fhat_cv=colSums(Fhat_cv0)/(nym-1)
+        
+        PRhatcv[w2,w1] = sum((dy_cv*(0.5 - Fhat_cv) > 0) +  0.5* ( dy_cv*(0.5 - Fhat_cv) ==0)) /nym
+        MSPEcv[w2,w1] =sum((diag(resid_cv)^2))/nym																			                                                                                
+      }
+      
+      
+    }
+    
+    
+    ##################################################################################################
+    # the optimal cluster number chosen by either maximizing PR(c) or minimizing MSPE(c) for c=      #
+    # 2,…,K                                                                                          #
+    ##################################################################################################
+    PR.MSPE = cbind(seq(1,(K-1)), MSPEcv[,l_m], PRhatcv[,l_m])
+    nclustera =  PR.MSPE[PR.MSPE[,2] == min(PR.MSPE[,2]),][1]  #nclustera + 1 = c
+    nclusterb =  PR.MSPE[PR.MSPE[,3] == max(PR.MSPE[,3]),][1]  #nclusterb + 1 = c
+    
+    
+    ##################################################################################################
+    # the information matrix for the stability of PR and MSPE for varying M_l, l = 6,…,m0            #
+    ##################################################################################################
+    rbind( PRhatcv[nclusterb, ] , MSPEcv[nclustera, ] )
+    
+    ##################################################################################################
+    # the LSEs of robot strengths in the proposed OPR model for varying Ml, l = 6,…,m0               #
+    ##################################################################################################
+    bhatinfm = matrix( 0 , K , l_m)
+    
+    for( w in 1:l_m){
+      ym = c( yr[-c( (c0[w] +1) : c0[l_m+1] )] , yb[-c( (c0[w] +1) : c0[l_m+1] )])
+      nym = length(ym)
+      dim(ym) = c(nym,1)
+      xm = rbind( xr[-c(c0[w] + 1 : c0[l_m+1]) ,] , xb[-c(c0[w] + 1 : c0[l_m+1]) ,]) 
+      
+      grp = cluster.inf0[(1:K),nclusterb]
+      grp0 = sort(unique(grp))
+      lgrp0 = length(grp0)
+      INF1 = outer( grp , grp0 , FUN="==") *1
+      bhatinfm[ ,w] = bhat.fn(xm,ym)
+    }
+    
+    
+    
+    # estimated strengths of robots in the playoff stage
+    bhatpo = cbind(ID,bhatinfm,bhat0,po)
+    bhatpo = bhatpo[bhatpo[,(l_m+4)] != 0 , 1:(l_m+3)]
+    potinfmb = bhatpo[,3:(l_m+3)]
+    
+    # ranks of robot strengths in decreasing order
+    bhatm0rk = bhatinfm[order(bhatinfm[,l_m] , decreasing = TRUE), 1:l_m]
+    
+    
+    
+    ##################################################################################################
+    # the information matrix for the stability of robot strengths ( all robots, top-8 robots, robots #
+    # in the playoff stage) for varying M_l, l = 6,…,m0                                              #                       
+    ##################################################################################################
+    RCs = matrix(0,3,l_m-1)
+    
+    for(w in 1:(l_m-1)){
+      RCs[1,w] = RC(bhatm0rk[,w] , bhatm0rk[,(w+1)])
+      RCs[2,w] = RC(bhatm0rk[(1:8),w] , bhatm0rk[(1:8),(w+1)])
+      RCs[3,w] = RC(potinfmb[,w] , potinfmb[,(w+1)])
+    }
+    
+    
+    ##################################################################################################
+    # PRs and MSPEs, which are computed based on playoff data with c-hat and K clusters              #  
+    ##################################################################################################
+    yhatpor = xrpo %*% bhatinfm[,l_m]
+    yhatpob = xbpo %*% bhatinfm[,l_m]
+    yhatpor0 = xrpo %*% bhat0
+    yhatpob0 = xbpo %*% bhat0
+    
+    Fhatpo =  colSums ( outer( (yr - yb) - (xr - xb) %*% bhatinfm[,l_m] , -( yhatpor - yhatpob ) , FUN = "<=")*1  ) / ( M)
+    dim(Fhatpo) = c(Mpo,1)
+    PRpo = (sum( ((yrpo - ybpo)*(0.5 - Fhatpo) > 0)*1 +  0.5* ((yrpo - ybpo)*(0.5 - Fhatpo) ==0)*1)) / (Mpo )
+    MSPEpo =  sum(((yrpo - ybpo) - (xrpo - xbpo) %*% bhatinfm[,l_m])^2)/ Mpo
+    
+    Fhatpo0 =  colSums ( outer( (yr - yb) - (xr - xb) %*% bhat0 , -( yhatpor0 - yhatpob0 ) , FUN = "<=")*1  ) / ( M)
+    dim(Fhatpo0) = c(Mpo,1)
+    PRpo0 = (sum( ((yrpo - ybpo)*(0.5 - Fhatpo0) > 0)*1 +  0.5* ((yrpo - ybpo)*(0.5 - Fhatpo0) ==0)*1)) / (Mpo )
+    MSPEpo0 =  sum(((yrpo - ybpo) - (xrpo - xbpo) %*% bhat0)^2)/ Mpo
+    
+    # information matrix for PRs and MSPEs in the qualification and playoff stages
+    rbind(c( PRhatcv[nclusterb,l_m] , PRhatcv[K-1,l_m] , MSPEcv[nclustera,l_m] , MSPEcv[ K-1 ,l_m] ) , c( PRpo , PRpo0 , MSPEpo , MSPEpo0) )
+
+    
+    ##################################################################################################
+    #agreement metrics (correlation, precision, and recall) between FRC ratings (all robots, top-8   #
+    #robots, robots in the playoff stage) and the proposed robots strength estimates                 #
+    ##################################################################################################
+    
+    # rank correlations computed under the OPR and OPRC models
+    RCm0 = matrix(0 , 2 , 3)
+    
+    RCm0[1,1] = RC( seq(K,1) , bhat.inf0[(1:K),l_n])
+    RCm0[1,2] = RC( seq(8,1) , bhat.inf0[(1:8),l_n])
+    RCm0[1,3] = RC( seq(Kpo,1) , bhatpo[,l_m+3])
+    RCm0[2,1] = RC( seq(K,1) , bhatinfm[,l_m])
+    RCm0[2,2] = RC( seq(8,1) , bhatinfm[1:8,l_m])
+    RCm0[2,3] = RC( seq(Kpo,1) , bhatpo[,l_m+2])
+    
+    # precisions and recalls computed under the OPR and OPRC models
+    Pm0 = order(bhat.inf0[(1:K),l_n],decreasing = TRUE)
+    Pm = order(bhatinfm[,l_m],decreasing = TRUE)
+    
+    Pre_Rec = matrix(0,K,4)
+    
+    for(w in 1:K){
+      Pre_Rec[w,1] = sum(outer(Pm0[1:w],seq(1,w),FUN="==")*1)/w
+      Pre_Rec[w,2] = sum(outer(Pm[1:w],seq(1,w),FUN="==")*1)/w
+      Pre_Rec[w,3] = sum(outer(Pm0[1:w],seq(1,8),FUN="==")*1)/8
+      Pre_Rec[w,4] = sum(outer(Pm[1:w],seq(1,8),FUN="==")*1)/8
+    }
+
+    
\ No newline at end of file