From b1d051caa65ef573ad094a55c7241a1d4e99981a Mon Sep 17 00:00:00 2001 From: 36cf6267d1f8856c85d87f0f287ea039 <36cf6267d1f8856c85d87f0f287ea039@app-learninglab.inria.fr> Date: Thu, 17 Jun 2021 16:12:52 +0000 Subject: [PATCH] Upload New File --- module4/Reproduction/OPRC1.R | 301 +++++++++++++++++++++++++++++++++++ 1 file changed, 301 insertions(+) create mode 100644 module4/Reproduction/OPRC1.R diff --git a/module4/Reproduction/OPRC1.R b/module4/Reproduction/OPRC1.R new file mode 100644 index 0000000..c8b8d1c --- /dev/null +++ b/module4/Reproduction/OPRC1.R @@ -0,0 +1,301 @@ +#--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- +# 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 -- 2.18.1