diff --git a/.DS_Store b/.DS_Store index 3e28009..bd63b96 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/DESCRIPTION b/DESCRIPTION index fe5d837..fb5fbf8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,11 +7,12 @@ Author@R: c(person("Ariel" "Shin", role-c("aut"), email="arishin@ucdavis.edu") ) Maintainer: Norm Matloff -Date: 12/12/2016 +Date: 5/27/2017 Title: Advanced Recommender System Description: Toolbox for recommender systems, including novel methods, - plotting, parallel computation, and so on. Requires Recosystem of at least 0.4 -Depends: partools, lme4, recosystem, shiny + plotting, parallel computation, and so on. +Depends: partools, lme4, shiny, dummies +Suggests: recosystem LazyLoad: no ByteCompile: yes NeedsCompilation: no diff --git a/NAMESPACE b/NAMESPACE index c85ea04..498428d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,9 +1,7 @@ S3method(predict,usrData) S3method(predict,RecoS3) - -importFrom("recosystem", Reco, data_memory, out_file) - +S3method(predict,ydotsMM) export(findYdotsMM,predict.ydotsMM,plot.ydotsMM,trainMM, findYdotsMLE,predict.ydotsMLE,trainMLE, @@ -12,7 +10,7 @@ export(findYdotsMM,predict.ydotsMM,plot.ydotsMM,trainMM, plot.ydotsMM,plot.xvalb, getTrainSet,getTestSet,trainReco,xvalReco, focusGrp, - trainMultiplic,predict.MMmultiplic,xvalMultiplic, - ratingness,covratingness) + ratingness,covratingness, + gassign,getInstEval,getML) diff --git a/R/.Rapp.history b/R/.Rapp.history deleted file mode 100644 index e69de29..0000000 diff --git a/R/.Rhistory b/R/.Rhistory deleted file mode 100644 index e69de29..0000000 diff --git a/R/Utils.R b/R/Utils.R index 108485e..7f89322 100644 --- a/R/Utils.R +++ b/R/Utils.R @@ -17,3 +17,68 @@ covratingness <- function(ratingsIn,splitCol) { tmp[ratingsIn[,splitCol]] } +# get the instructor evaluation data, and set it up, including +# covariates + +getInstEval <- function() { + data(InstEval) + ivl <- InstEval + # convert from factors + ivl$s <- as.numeric(ivl$s) + ivl$d <- as.numeric(ivl$d) + ivl$studage <- as.numeric(ivl$studage) + ivl$lectage <- as.numeric(ivl$lectage) + ivl$service <- as.numeric(ivl$service) + # make correct format, choose + ivl <- ivl[,c(1,2,7,3:6)] + # create dummy variables in place of dept + library(dummies) + dms <- dummy(ivl$dept) + dms <- as.data.frame(dms) + dms$dept2 <- NULL + ivl$dept <- NULL + ivl <- cbind(ivl,dms) + gassign('ivl','ivl') +} + +# get the MovieLense 100K evaluation data, and set it up, including +# covariates + +# if needDownload, then download will be done; datadir is the directory +# containing the data, assumed by default to be in the current working +# directory + +getML <- function(needDownload=FALSE,datadir='./ml-100k') { + if (needDownload) { + # 5 Mb + download.file( + 'http://files.grouplens.org/datasets/movielens/ml-100k.zip', + 'ml-100k.zip') + unzip('ml-100k.zip') + } + currdir <- getwd() # leave a trail of bread crumbs + setwd(datadir) + ud <- read.table('u.data',header=F,sep='\t') + uu <- read.table('u.user',header=F,sep='|') + ui <- read.table('u.item',header=F,sep='|') + setwd(currdir) # follow the trail back + ud <- ud[,-4] # remove timestamp, leaving user, item, rating + uu <- uu[,1:3] # user, age, gender + ui <- ui[,c(1,6:24)] # item num, genres + names(ud) <- c('user','item','rating') + names(uu) <- c('user','age','gender') + names(ui)[1] <- 'item' + names(ui)[-1] <- gsub('V','GNR',names(ui)[-1]) # GNR = genre + uu$gender <- as.integer(uu$gender == 'M') + uduu <- merge(ud,uu) + uduuui <- merge(uduu,ui) + gassign('uduu','uduu') + gassign('uduuui','uduuui') +} + +# not-so-stealth global assignment +gassign <- function(gname,rhsname) { + cmd <- paste(gname,'<<-',rhsname) + eval(parse(text=cmd),envir=parent.frame()) +} + diff --git a/R/cosine.R b/R/cosine.R index c918b37..aeb4e0b 100644 --- a/R/cosine.R +++ b/R/cosine.R @@ -35,6 +35,7 @@ predict.usrData <- function(origData,newData,newItem, # (NA,NA); defined for use by sapply() below checkNewItem <- function(oneUsr) { tmp <- match(oneUsr$itms, newItem) + #message(oneUsr) if (all(is.na(tmp))) { c(NA,NA) } diff --git a/R/findUsrItmData.R b/R/findUsrItmData.R index d5402eb..11911d5 100644 --- a/R/findUsrItmData.R +++ b/R/findUsrItmData.R @@ -78,13 +78,13 @@ formUserData <- function(ratingsIn,usrCovs=NULL,itmCats=NULL,fileOut='') { retval } -### # construct a new object of class 'usrDatum' -### -### formUserDatum <- function(itms,ratings,userID=NULL) { -### obj <- list(itms = itms, ratings=ratings,userID=userID) -### class(obj) <- 'usrDatum' -### obj -### } +# construct a new object of class 'usrDatum' + +formUserDatum <- function(itms,ratings,userID=NULL) { + obj <- list(itms = itms, ratings=ratings,userID=userID) + class(obj) <- 'usrDatum' + obj +} # utility: find input row for a given user, item findInputRow <- function(ratingsIn,usrID,itmID) { diff --git a/R/findYdotsMM.R b/R/findYdotsMM.R index abe38f6..1bb4aba 100644 --- a/R/findYdotsMM.R +++ b/R/findYdotsMM.R @@ -1,57 +1,62 @@ +# modified version of earlier findYdotsMM() etc., May 27, 2017, with +# new approach to use of covariates U_ijk for user i, item j; +# NOTE: the covariates will be centered + +# basic model: +# +# Y_ij = +# mu + sum_k gamma_k U_ijk + alpha_i + beta_j + eps +# +# that means + +# E(Y_ij | U_ik, V_jk) = sum_k gamma_k U_ik + sum_k delta_k V_jk + +# and the coefficients can be estimated via lm() without a const term + # arguments: -# ratingsIn: input data, with first cols (userID,itemID,rating, -# covariates); data frame, unless cls is non-null, in which -# case this argument is the quoted name of the distributed -# data frame -# regressYdots; if TRUE, apply lm() to the estimated latent factors -# (and their product), enabling rating prediction from the -# resulting linear function of the factors; currently -# only implemented if have no covariates -# cls: an R 'parallel' cluster +# ratingsIn: input data, with cols (userID,itemID,rating, +# covariates); data framee # value: # S3 class of type "ydots", with components: -# Y..: grand mean (0 if have covariates) -# Yi.: vector of mean ratings for each user -# Y.j: vector of mean ratings for each item -# regOjb: if have covariates, regression output, e.g. coefs - -findYdotsMM <- function(ratingsIn,regressYdots=FALSE,cls=NULL) { - require(partools) - if (is.null(cls)) { - users = ratingsIn[,1] - items = ratingsIn[,2] - ratings = ratingsIn[,3] - nms <- names(ratingsIn) - haveCovs = ncol(ratingsIn) > 3 - if (haveCovs) { - frml = as.formula(paste(nms[3],'~ .')) - lmout = lm(frml,data=ratingsIn[,-(1:2)]) - fits = lmout$fitted.values - ratings = ratings - fits - Y.. = 0 - } else Y.. = mean(ratings) - Yi. = tapply(ratings,users,mean) # means of all ratings per user - Y.j = tapply(ratings,items,mean) # means of all ratings per item - } else { - stop('parallel version under construction') - } - ydots = list(grandMean=Y..,usrMeans=Yi.,itmMeans=Y.j) +# Y..: grand mean, est of mu +# Yi.: vector of mean ratings for each user, ests. of alpha_i +# Y.j: vector of mean ratings for each item, ests. of betaa_j +# regObj: if have covariates, regression output, e.g. coefs + +findYdotsMM <- function(ratingsIn,regressYdots=FALSE) { + users <- ratingsIn[,1] + items <- ratingsIn[,2] + # IMPORTANT NOTE: + # user and item IDs may not be consecutive; even if they are + # consecutive in the original, if we do cross-validation, this + # may not be the case; so switch to character IDs + users <- as.character(users) + items <- as.character(items) + ratings <- ratingsIn[,3] + nms <- names(ratingsIn) + Y.. <- mean(ratings) + Yi. <- tapply(ratings,users,mean) # means of all ratings per user + Y.j <- tapply(ratings,items,mean) # means of all ratings per item + ydots <- list(grandMean=Y..,usrMeans=Yi.,itmMeans=Y.j) + haveCovs <- ncol(ratingsIn) > 3 if (haveCovs) { - ydots$regObj = lmout + # center the covs + tmp <- scale(ratingsIn[,-(1:3)],scale=FALSE) + ratingsIn[,-(1:3)] <- tmp + ydots$covmeans <- attr(tmp,'scaled:center') + # regress, no constant term; could do a weighted least squares, + # using the Ni, but since the latter is random too, not needed + frml <- as.formula(paste(nms[3],'~ .-1')) + lmout <- lm(frml,data=ratingsIn[,-(1:2)]) + ydots$lmout <- lmout + Ni <- tapply(ratings,users,length) # number of ratings per user + ydots$Ni <- Ni } - if (regressYdots && !haveCovs) { - yi. = Yi.[ratingsIn[,1]] - y.j = Y.j[ratingsIn[,2]] - # yij = yi. * y.j - # ydots$regressYdots = coef(lm(ratings ~ yi. + y.j + yij)) - ydots$regressYdots = coef(lm(ratings ~ yi. + y.j)) - } - class(ydots) = 'ydotsMM' invisible(ydots) } @@ -59,33 +64,152 @@ findYdotsMM <- function(ratingsIn,regressYdots=FALSE,cls=NULL) { # alias trainMM <- findYdotsMM -# check +# predict() method for the 'ydotsMM' class + +# in predicting for user i, the code looks at N_i, the number of ratings +# by user i; if that number is below minN, the prediction comes from +# user i's covariate information instead of from Yi. + +# arguments + +# testSet: data frame in same form as ratingsIn above except that there +# is no ratings column; thus covariates, if any, are shifted +# leftward one slot, i.e. userID, itemID, cov1, cov2... +# ydotsObj: the output of findYdotsMM() +# minN: if Ni < minN and have covariates, use the latter instead of Yi; +# see above + +# returns vector of predicted values for testSet +predict.ydotsMM = function(ydotsObj,testSet,minN=0) { + haveCovs <- ncol(testSet) > 2 + # see comment on as.character() above + ts1 <- as.character(testSet[,1]) # user IDs, char form + # tmp will basically consist of the user means, except that in the + # covariate case some will be replaced by predict.lm() values + Y.. + if (!haveCovs) { + tmp <- ydotsObj$usrMeans[ts1] + } + else { + colmeans <- ydotsObj$covmeans + testSet[,-(1:2)] <- + scale(testSet[,-(1:2)],center=colmeans,scale=FALSE) + tmp <- + ifelse (ydotsObj$Ni[ts1] >= minN, + ydotsObj$usrMeans[ts1], + predict(ydotsObj$lmout,testSet[,-(1:2),drop=FALSE]) + + ydotsObj$grandMean + ) + } + testSet$pred <- tmp + + ydotsObj$itmMeans[as.character(testSet[,2])] - ydotsObj$grandMean + testSet$pred +} + +# test case checkyd <- function() { - check <- - data.frame(userID = c(1,3,2,1,2),itemID = c(1,1,3,2,3),ratings=6:10) + check <- data.frame( + userID <- c(1,3,2,1,2), + itemID <- c(1,1,3,2,3), + ratings <- 6:10) + names(check) <- c('u','i','r') print(check) print(findYdotsMM(check)) check$cv <- c(1,4,6,2,10) + names(check)[4] <- 'x' print(check) - print(findYdotsMM(check)) + cout <- findYdotsMM(check) + print(cout) + testset <- check[1:2,-3] + testset$x <- c(5,8) + print(predict(cout,testset,2)) } -# predict from output from coef(lm()); not used currently -predict.lmcoef = function(coefs,testSet) { - cbind(1,as.matrix(testSet)) %*% as.vector(coefs) + + +predict.NM = function(ydotsObj,testSet) { + testSet$pred = ydotsObj$usrMeans[as.character(testSet[,1])] + + ydotsObj$itmMeans[as.character(testSet[,2])] - ydotsObj$grandMean +# if (!is.null(ydotsObj$regObj)) + # testSet$pred = testSet$pred + predict(ydotsObj$regObj,testSet[,-(1:2)]) + testSet$pred } -# predict() method for the 'ydots' class -# -# testSet in same form as ratingsIn in findYdots(), except that there -# is no ratings column; regObj is as in the output of findYdots() -# -# returns vector of predicted values for testSet -predict.ydotsMM = function(ydotsObj,testSet) { - testSet$pred = ydotsObj$usrMeans[as.character(testSet[,1])] + - ydotsObj$itmMeans[as.character(testSet[,2])] - ydotsObj$grandMean - if (!is.null(ydotsObj$regObj)) - testSet$pred = testSet$pred + predict(ydotsObj$regObj,testSet[,-(1:2)]) - testSet$pred +getDiff <- function(fullMatrix, filledNMF){ + indicies <- which(!is.na(fullMatrix), arr.ind = TRUE) + diff <- matrix(data = NA,nrow = nrow(indicies), 1) + for(i in 1:nrow(indicies)) + for(j in 1:ncol(indicies)){ + diff[i,1] <- round(fullMatrix[indicies[i,1],indicies[j,2]]) - filledNMF[indicies[i,1],indicies[j,2]] + } + + diff +} + + +getAcc <- function(fullMatrix,filledNMF, threshold = 0.5) { + diff <- abs(getDiff(fullMatrix, filledNMF)) + avgDiff <- colMeans(abs(diff), na.rm = TRUE, dims =1) # average difference is 1.19 + # Remove the na's from diff + # Calculate the values less than the threshold + count = 0; + tester <- na.exclude(diff) + for(i in 1:nrow(diff)) + { + if(diff[i] < 0.5 ) + count = count + 1 + } + acc <- count/ nrow(diff) + acc + +} + +train.NM <- function(ratingsIn, trainprop = 0.5,cls = NULL, + rnk = 10, recosystem = FALSE,regressYdots=FALSE) + { + require(NMF) + library(NMF) + if(recosystem == TRUE){ + source((paste(system.file(package = 'rectools'), + 'recosys/xvalRecoParallel.R', sep = ""))) + res <- xvalReco(ratingsIn,trainprop,cls,rnk) + } + else { + origMatrix <- buildMatrix(ratingsIn) # Matrix A (Step 1) + fullMatrix <- origMatrix + approxMatrix <- findYdotsMM(ratingsIn) # Matrix V (Step 2) + + naMatrix <- as.data.frame(which(is.na(fullMatrix) == TRUE, arr.ind = TRUE)) + naMatrix$ratings <- NA + + preds <- predict.NM(approxMatrix,naMatrix) # Step 3 + + fullMatrix[which(is.na(fullMatrix))] <- preds + fullMatrix[fullMatrix < 0] <- 0 + require(NMF) + result <- nmf(as.matrix(fullMatrix),10) # Step 4 + filledNMF <- as.matrix(result@fit@W) %*% as.matrix(result@fit@H) + avgd <- getAcc(fullMatrix, filledNMF) # 1.149 for the movielens 100k data + print(avgd) + + } + filledNMF +} + + +buildMatrix <- function(ratingsIn,NAval=NA){ + # deal with possible factors + dmax <- function(d) { + if (is.factor(d)) return(length(levels(d))) + max(d) + } + users = ratingsIn[,1] + movies = ratingsIn[,2] + rating = ratingsIn[,3] + newMatrix = matrix(NAval, + nrow = dmax(users), ncol = dmax(movies)) + for(rows in 1:nrow(ratingsIn)){ + newMatrix[ratingsIn[rows,1],ratingsIn[rows,2]] = ratingsIn[rows,3] + } + return (newMatrix) } diff --git a/R/findYdotsMM.old b/R/findYdotsMM.old new file mode 100644 index 0000000..38f4f30 --- /dev/null +++ b/R/findYdotsMM.old @@ -0,0 +1,143 @@ + +# arguments: + +# ratingsIn: input data, with first cols (userID,itemID,rating, +# covariates); data frame, unless cls is non-null, in which +# case this argument is the quoted name of the distributed +# data frame +# regressYdots; if TRUE, apply lm() to the estimated latent factors +# (and their product), enabling rating prediction from the +# resulting linear function of the factors; currently +# only implemented if have no covariates +# cls: an R 'parallel' cluster + +# value: + +# S3 class of type "ydots", with components: + +# Y..: grand mean (0 if have covariates) +# Yi.: vector of mean ratings for each user +# Y.j: vector of mean ratings for each item +# regOjb: if have covariates, regression output, e.g. coefs + +findYdotsMM <- function(ratingsIn,regressYdots=FALSE,cls=NULL) { + require(partools) + if (is.null(cls)) { + users = ratingsIn[,1] + items = ratingsIn[,2] + ratings = ratingsIn[,3] + nms <- names(ratingsIn) + haveCovs = ncol(ratingsIn) > 3 + if (haveCovs) { + frml = as.formula(paste(nms[3],'~ .')) + lmout = lm(frml,data=ratingsIn[,-(1:2)]) + fits = lmout$fitted.values + ratings = ratings - fits + Y.. = 0 + } else Y.. = mean(ratings) + Yi. = tapply(ratings,users,mean) # means of all ratings per user + Y.j = tapply(ratings,items,mean) # means of all ratings per item + } else { + stop('parallel version under construction') + } + ydots = list(grandMean=Y..,usrMeans=Yi.,itmMeans=Y.j) + if (haveCovs) { + ydots$regObj = lmout + } + if (regressYdots && !haveCovs) { + yi. = Yi.[ratingsIn[,1]] + y.j = Y.j[ratingsIn[,2]] + # yij = yi. * y.j + # ydots$regressYdots = coef(lm(ratings ~ yi. + y.j + yij)) + ydots$regressYdots = coef(lm(ratings ~ yi. + y.j)) + } + + class(ydots) = 'ydotsMM' + invisible(ydots) +} + +# alias +trainMM <- findYdotsMM + +# check +checkyd <- function() { + check <- + data.frame(userID = c(1,3,2,1,2),itemID = c(1,1,3,2,3),ratings=6:10) + print(check) + print(findYdotsMM(check)) + check$cv <- c(1,4,6,2,10) + print(check) + print(findYdotsMM(check)) +} + +# predict from output from coef(lm()); not used currently +predict.lmcoef = function(coefs,testSet) { + cbind(1,as.matrix(testSet)) %*% as.vector(coefs) +} + +# predict() method for the 'ydots' class +# +# testSet in same form as ratingsIn in findYdots(), except that there +# is no ratings column; regObj is as in the output of findYdots() +# +# returns vector of predicted values for testSet +predict.ydotsMM = function(ydotsObj,testSet) { + testSet$pred = ydotsObj$usrMeans[as.character(testSet[,1])] + + ydotsObj$itmMeans[as.character(testSet[,2])] - ydotsObj$grandMean + if (!is.null(ydotsObj$regObj)) + testSet$pred = testSet$pred + predict(ydotsObj$regObj,testSet[,-(1:2)]) + testSet$pred +} +# buildMatrix(): A function that changes ratingsIn into a mostly empty matrix upon which we fill +buildMatrix <- function(ratingsIn,NAval=NA){ + # deal with possible factors + dmax <- function(d) { + if (is.factor(d)) return(length(levels(d))) + max(d) + } + users = ratingsIn[,1] + movies = ratingsIn[,2] + rating = ratingsIn[,3] + newMatrix = matrix(NA, + nrow = dmax(users), ncol = dmax(movies)) + for(rows in 1:nrow(ratingsIn)){ + newMatrix[ratingsIn[rows,1],ratingsIn[rows,2]] = ratingsIn[rows,3] + } + return (newMatrix) +} + +# Wrapper for recosystem +# TrainNM will check if recosystem is called, if so, call recosystem function +# If not, build ratingsIn into a mostly matrix which would be users x items dimensions +# Approximate the unknown ratings using findYdotsMM +# Substitute those values and call NMF() +trainNM <- function(ratingsIn, trainprop = 0.5,cls = NULL, + rnk = 10, recosystem = FALSE,regressYdots=FALSE){ + require(NMF) + library(NMF) + if(recosystem == TRUE){ + source((paste(system.file(package = 'rectools'), + 'recosys/xvalRecoParallel.R', sep = ""))) + res <- xvalReco(ratingsIn,trainprop,cls,rnk) + } + else { + fullMatrix <- buildMatrix(ratingsIn) # Matrix A (Step 1) + fullMatrix[which(fullMatrix == 0)] = NA + approxMatrix <- findYdotsMM(ratingsIn) # Matrix V (Step 2) + + naMatrix <- as.data.frame(which(is.na(fullMatrix) == TRUE, arr.ind = TRUE)) + naMatrix$ratings <- NA + + + preds <- predict.ydotsMM(approxMatrix,naMatrix) # Step 3 + + + + fullMatrix[which(is.na(fullMatrix))] <- preds + fullMatrix[fullMatrix < 0] <- 0 + require(NMF) + res <- nmf(as.matrix(fullMatrix),10) # Step 4 + } + res +} + diff --git a/R/hybrid.R b/R/hybrid.R new file mode 100644 index 0000000..79a215a --- /dev/null +++ b/R/hybrid.R @@ -0,0 +1,129 @@ +#predict using the k nearest neighbour from the "basis" of the user-space(origData) +# covariates (e.g. age, gender) and item type preferences (e.g. +# preferred movie genres) are allowed + +########################## predict() function ########################## + +# predict newData argument from origData argument using the hybrid kNN-Matrix Product Method + +# arguments: + +# datain: training set, object of class 'usrData' (see file +# findUsrItmData.R) +# newData: data point (just one for now) to be predicted, object of +# class 'usrDatum' +# newItem: ID of the item rating to be predicted for the user in +# newData +# wtcovs: weight to put on covariates; NULL if no covs +# wtcats: weight to put on item categories; NULL if no cats +# k: a list of the number of nearest neigbhors to be considered while predicting +# r: rank for factorisation + +# value: + +# predicted ratings for newData + +hybdpredict.usrData <- function(datain,newData,newItem, + k,r, wtcovs=NULL,wtcats=NULL) { + + origData <- formUserData(datain) + #getting the "basis" matix of the "user space" + if (r >= min(dim(origData))) { + stop('no data on this item') + } + else{ + appMatrix = Reco() + train_set = data_memory(unlist(datain[1]), unlist(datain[2]), + unlist(datain[3]), index1 = TRUE) + appMatrix$train(train_set) + #get factors after NMF -- need to use $tune() fix later + factors <- appMatrix$output(out_memory(), out_memory()) + #P- is the "user space basis"-- see RecoSys for P and Q + P <- factors$P + usrs <- rep(1:dim(P)[1], each = dim(P)[2]) + itms <- rep.int(1:dim(P)[2], dim(P)[1]) + #get ratings from "user space basis" + rtngs <- as.vector(t(P)) + #form basis (is of type usrData) + basis <- formUserData(data.frame(usrs, itms, rtngs)) + } + + # action of checkNewItem(): here oneUsr is one representative-user record in + # the "basis" of the use-space; + #the function will check which item in the record, if any, + # equals newItem, and will return that value and the corresponding + # rating; defined for use by sapply() below + checkNewItem <- function(oneUsr) { + #itms <- as.list(as.data.frame(t(oneUsr$itms))) + tmp <- match(oneUsr$itms[[1]], newItem) + if (all(is.na(tmp))) { + c(NA, NA) + } + else + { + whichOne <- which(!is.na(tmp)) + c(whichOne,oneUsr$ratings[whichOne]) + } + } + found <- as.matrix(sapply(basis,checkNewItem)) + # found is of dimensions 2 x number of users; + # found[1,i] = j means origData[[i]]$itms[j] = newItem; + # found[1,i] = NA means newItem wasn't rated by user i;ß + # found[2,i] will be the rating in the non-NA case + # we need to get rid of the NA users + whoHasIt <- which(!is.na(found[1,])) + # whoHasIt[i] is the index, i.e. user ID, of the i-th user who has + # rated newData + if (is.null(whoHasIt) | length(whoHasIt)==0) return(NA) # no one rated this item + basis <- basis[whoHasIt] + # now basis only has the relevant users, the ones who have rated + # newItem, so select only those columns of the found matrix + found <- found[,whoHasIt,drop=FALSE] + #found <- found[!is.na(found)] + # find the distance from newData to one user y of origData; defined for + # use in sapply() below + onecos <- function(y) cosDist(newData,y,wtcovs,wtcats) + cosines <- sapply(basis,onecos) + # the vector cosines contains the distances from newData to all the + # original data points; + + #action of findKnghbourRtng(): predict rating based on each k[i] neighbours + #x = k[i] + #if x > neighbours present in the dataset, then the maximum number of meighbours is used + findKnghbourRtng <- function(x){ + #x can be atmost the number of neighbours in the dataset + x <- min(x, length(cosines)) + #klarge is a vector containing the indices of the x closest neighbours + klarge <- order(cosines,decreasing=TRUE)[1:x] + mean(as.numeric(found[2, klarge])) + } + sapply(k, findKnghbourRtng) +} + +# cosDist() find cosine distance between x and y, elements of an object +# of 'usrData' class; only items rated in both x and y are used; if none +# exist, then return NaN +cosDist <- function(x,y,wtcovs,wtcats) { + # rated items in common + commItms <- intersect(x$itms[[1]],y$itms[[1]]) + if (is.null(commItms)| length(commItms)==0) return(NaN) + # where are they in x and y? + xwhere <- which(!is.na(match(x$itms[[1]],commItms))) + ywhere <- which(!is.na(match(y$itms[[1]],commItms))) + xrats <- x$ratings[xwhere,] + yrats <- y$ratings[ywhere,] + cosTot <- xrats[[1]] %*% yrats[[1]] + xl2 <- sum(xrats^2) + yl2 <- sum(yrats^2) + if (!is.null(wtcovs)) { + cosTot <- cosTot + wtcovs * x$cvrs %*% y$cvrs + xl2 <- xl2 + sum((wtcovs*x$cvrs)^2) + yl2 <- yl2 + sum((wtcovs*y$cvrs)^2) + } + if (!is.null(wtcats)) { + cosTot <- cosTot + wtcats * x$cats %*% t(y$cats) + xl2 <- xl2 + sum((wtcovs*x$cats)^2) + yl2 <- yl2 + sum((wtcovs*y$cats)^2) + } + cosTot / sqrt(xl2 * yl2) +} diff --git a/R/xvalMM.R b/R/xvalMM.R index 6bb0beb..8e83385 100644 --- a/R/xvalMM.R +++ b/R/xvalMM.R @@ -16,16 +16,15 @@ # (and their product), enabling rating prediction from the # resulting linear function of the factors; currently # only implemented if have no covariates -# cls: if non-null, +# minN: see findYdotsMM.R # value: # accuracy value xvalMM <- function(ratingsIn, trainprop=0.5, - regressYdots=FALSE,cls=NULL){ - if(!is.null(cls)) stop('parallel version under construction') - if(is.null(cls)) ratIn = ratingsIn else ratIn = get(ratingsIn) + regressYdots=FALSE,minN=0){ + ratIn = ratingsIn # split into random training and validation sets nrowRatIn = nrow(ratIn) rowNum = floor(trainprop * nrowRatIn) @@ -35,13 +34,13 @@ xvalMM <- function(ratingsIn, trainprop=0.5, trainItems = trainingSet[,2] trainUsers = trainingSet[,1] # get means - means = findYdotsMM(trainingSet,regressYdots,cls) + means = findYdotsMM(trainingSet,regressYdots) # Y.. = means$grandMean # Yi. = means$usrMeans # Y.j = means$itmMeans testIdxs = setdiff(1:nrowRatIn,trainIdxs) testA = ratIn[testIdxs,] - testA$pred = predict(means,testA[,-3]) # predict.ydotsMM + testA$pred = predict(means,testA[,-3],minN) # predict.ydotsMM # need to integrate the following into predict.ydotsMM if (regressYdots) { yi. = means$usrMeans[testA[,1]] diff --git a/README.md b/README.md index ba45ef3..a106b29 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,8 @@ * Parallel computation. -* Novel variations on statistical latent factor model. +* Novel variations on common models, e.g. a hybrid of + NMF and k-Nearest Neighbor. * Plotting. @@ -27,9 +28,9 @@ A simple random effects latent factor model is E(Y) = μ + αi + βj -where Yij is the rating, α and βj being -specific latent effects for user i and item j, e.g. movie reviewer i and -movie j. +where Yij is the rating, with αi and +βj being specific latent effects for user i and item j, +e.g. movie reviewer i and movie j. Though typically Maximum Likelihood Estimation is used for latent factor models, this is computationally infeasible on large data sets. Instead, @@ -41,22 +42,23 @@ overall mean Y.. The predicted value of Yij is then Yi. + Y.j - Y.. -Computation is simple, conducted by our function **findYdotsMM**; -prediction is done on the output by our function **predict.ydotsMM**. +Computation is simple, with estimation conducted by our function +**findYdotsMM()**; prediction is done on the output by our function +**predict.ydotsMM()**. A novel enhancement in the package is to allow for different weights to be given to the αi and βj components in -the MM version. +the MM version. (With MLE it wouldn't matter, just changing the +variances.) We do make MLE available. Here αi and βj are assumed to have independent normal distributions -with different variances. (The error term -εij = Yij - EYij -is assumed independent of αi and βj, -with variance constant across i and j.) -We piggyback R's **lme4** package, -forming a wrapper for our application, and adding our function -**predict.ydotsMLE**, also a wrapper suited for our context. Since MLE +with different variances. (The error term εij = +Yij - EYij is assumed independent of +αi and βj, with variance constant +across i and j.) We piggyback R's **lme4** package, forming a wrapper +for our application, and adding our function **predict.ydotsMLE()** for +prediction, also an **lme4** wrapper suited for our context. Since MLE computation can be voluminous, our package offers a parallel version. Covariates are allowed for both the MM and MLE versions. @@ -75,7 +77,7 @@ Covariates are allowed for both the MM and MLE versions. > ivl$studage <- as.numeric(ivl$studage) > ivl$lectage <- as.numeric(ivl$lectage) > ivl$service <- as.numeric(ivl$service) -> # Make correct format, choose covs: +> # Make correct format (user ID, item ID, rating), choose covs: > ivl <- ivl[,c(1,2,7,3:6)] > # Create dummy variables in place of dept: > library(dummies) @@ -107,7 +109,7 @@ Covariates are allowed for both the MM and MLE versions. Let A denote the matrix of ratings, with Yij in row i, column j. Most of A is unknown, and we wish to predict the unknown values. -*Nonnegative matrix factorization* does this as follows: +Nonnegative Matrix Factorization (NMF) does this as follows: We find nonnegative matrices W and H, each of rank k, such that A is approximately equal to the product WH. Here k is a user-defined tuning diff --git a/inst/.DS_Store b/inst/.DS_Store index bb989f7..c1e0b87 100644 Binary files a/inst/.DS_Store and b/inst/.DS_Store differ diff --git a/R/findAlphBet.R b/inst/findAlphBet.R similarity index 93% rename from R/findAlphBet.R rename to inst/findAlphBet.R index 3175b37..744c74f 100644 --- a/R/findAlphBet.R +++ b/inst/findAlphBet.R @@ -2,9 +2,10 @@ # multiplicative random effects model: the rating Y, 0 or 1, is # generated by # -# EY_ij = q alpha_i beta_j +# EY_ij = b0 = b1 alpha_i beta_j -# for unobservable alpha, beta, and an unknown constant q +# for unobservable alpha, beta having means 1.0, +# and unknown constants b0, b1 # this is a very simple model, whose virtue is simplicity and quick # computation diff --git a/R/findAlphBet.R.iter b/inst/findAlphBet.R.iter similarity index 100% rename from R/findAlphBet.R.iter rename to inst/findAlphBet.R.iter diff --git a/inst/rectoolShiny2/.DS_Store b/inst/rectoolShiny2/.DS_Store index 5008ddf..1a71cb0 100755 Binary files a/inst/rectoolShiny2/.DS_Store and b/inst/rectoolShiny2/.DS_Store differ diff --git a/vignettes/rectools.Rnw b/inst/vignettes/rectools.Rnw similarity index 100% rename from vignettes/rectools.Rnw rename to inst/vignettes/rectools.Rnw diff --git a/R/xvalMultiplic.R b/inst/xvalMultiplic.R similarity index 100% rename from R/xvalMultiplic.R rename to inst/xvalMultiplic.R diff --git a/man/findYdots.Rd b/man/findYdots.Rd index 33430f4..6f94826 100644 --- a/man/findYdots.Rd +++ b/man/findYdots.Rd @@ -16,7 +16,7 @@ capability. Various plots can be displayed. } \usage{ -findYdotsMM(ratingsIn,regressYdots=FALSE,cls=NULL) +findYdotsMM(ratingsIn,regressYdots=FALSE,minN) trainMM(ratingsIn,regressYdots=FALSE,cls=NULL) predict.ydotsMM(ydotsObj,testSet) plot.ydotsMM(ydotsObj,ratingsIn) @@ -31,6 +31,9 @@ predict.ydotsMLE(ydotsObj,testSet) \item{regressYdots}{If TRUE, apply linear regression to the latent factors.} \item{cls}{An R \code{parallel} cluster.} + \item{minN}{If a prediction is to be made on a user with fewer than + this number of ratings and if there are covariates, predict from + the covariates.} \item{ydotsObj}{An object of class \code{'ydotsMM'} or \code{'ydotsMLE'}.} \item{testSet}{A data frame consisting of cases to be predicted. Format is the same as \code{ratingsIn}, except that there is no @@ -110,6 +113,11 @@ ivl <- ivl[,c(1,2,7,3:6)] library(dummies) dms <- dummy(ivl$dept) dms <- as.data.frame(dms) +# numeric names won't work, so change +dnms <- names(dms) +for (i in 1:length(dnms)) + dnms[i] <- paste('x',dnms[i],sep='') +names(dms) <- dnms dms$dept2 <- NULL ivl$dept <- NULL ivl <- cbind(ivl,dms) diff --git a/man/xvalReco.Rd b/man/xvalReco.Rd index b09d4f3..e6409e8 100644 --- a/man/xvalReco.Rd +++ b/man/xvalReco.Rd @@ -16,7 +16,7 @@ computation capability in some cases. } \usage{ - xvalMM(ratingsIn, trainprop = 0.5, regressYdots = FALSE, cls = NULL) + xvalMM(ratingsIn, trainprop = 0.5, regressYdots = FALSE, minN = 0) xvalMLE(ratingsIn, trainprop = 0.5, cls = NULL) xvalReco(ratingsIn, trainprop = 0.5,cls = NULL,rnk = 10) xvalCos(ratingsIn, k, usrCovs = NULL, itmCats = NULL, wtcovs = NULL, @@ -30,6 +30,10 @@ computation capability in some cases. \arguments{ \item{ratingsIn}{Input data frame. Within-row format is UserID, ItemID, rating} + \item{minN}{Applies to situations in which covariates are + present. In predicting for user i, then either Yi. or the + regression-based prediction will be used, depending on whether + N_i >= minN.} \item{trainprop}{The fraction of ratingsIn we want to use for our training set} \item{cls}{R \code{paralell} cluster.}