Data Mining Tutorial

Welcome to the companion website of  the “Mining Big Data to Extract Patterns and Predict Real-Life Outcomes” article published in Psychological Methods in 2016. 

Read In and Preprocess the Data

First, you need to download the sample data set. Next, unzip it. Then, you can read it into R: 

# Read in files (provide a full path, e.g., "~/Desktop/sample_dataset/users-likes.csv")

users <- read.csv("users.csv")

likes <- read.csv("likes.csv")

ul <- read.csv("users-likes.csv")    

# You can check what's inside each object using the following set of commands:

head(users)

head(likes)

head(ul)

tail(ul)

tail(users)

tail(likes)

dim(ul)

dim(users)

dim(likes)

# Match entries in ul with users and likes dictionaries

ul$user_row<-match(ul$userid,users$userid)

ul$like_row<-match(ul$likeid,likes$likeid)

# and inspect what happened: 

head(ul)

Next, let's build user-footprint Matrix M:

# Install Matrix library - run only once

install.packages("Matrix")

# Load Matrix library

require(Matrix)

# Construct the sparse User-Like Matrix M

M <- sparseMatrix(i = ul$user_row, j = ul$like_row, x = 1)

# Check the dimensions of M

dim(M)

# Save user IDs as row names in M

rownames(M) <- users$userid

# Save Like names as column names in M

colnames(M) <- likes$name

# Remove ul and likes objects (they won't be needed)

rm(ul, likes)

Next, trim matrix M to remove users and Likes that are too rare.

# Remove users/Likes occurring less than 50/150 times

repeat {                                       # Repeat whatever is in the brackets

  i <- sum(dim(M))                             # Check the size of M

  M <- M[rowSums(M) >= 50, colSums(M) >= 150]  # Retain only these rows/columns that meet the threshold

  if (sum(dim(M)) == i) break                  # If the size has not changed, break the loop

  }

# Check the new size of M

dim(M)

# Remove the users from users object that were removed

# from M

users <- users[match(rownames(M), users$userid), ]

# Check the new size of users

dim(users)

Reduce the Dimensionality of M

Using SVD:

# Preset the random number generator in R 

# for the comparability of the results

set.seed(seed = 68)

# Install irlba package (run only once)

install.packages("irlba")

# Load irlba and extract 5 SVD dimensions

library(irlba)

Msvd <- irlba(M, nv = 5)

# User SVD scores are here:

u <- Msvd$u

# Like SVD scores are here:

v <- Msvd$v

# The scree plot of singular values:

plot(Msvd$d)

You should also varimax-rotate the resulting SVD space.

# First obtain rotated V matrix:

# (unclass function has to be used to save it as an 

# object of type matrix and not loadings)

v_rot <- unclass(varimax(Msvd$v)$loadings)

# The cross-product of M and v_rot gives u_rot:

u_rot <- as.matrix(M %*% v_rot)

Using LDA:

# Install topicmodels package (run only once)

install.packages("topicmodels")

# Load it

library(topicmodels)

# Conduct LDA analysis, see text for details on setting

# alpha and delta parameters. 

# WARNING: this may take quite some time!

Mlda <- LDA(M, control = list(alpha = 10, delta = .1, seed=68), k = 5, method = "Gibbs")

# Extract user LDA cluster memberships

gamma <- Mlda@gamma

# Extract Like LDA clusters memberships

# betas are stored as logarithms, 

# function exp() is used to convert logs to probabilities

beta <- exp(Mlda@beta) 

# Log-likelihood of the model is stored here:

Mlda@loglikelihood

# and can be also accessed using logLik() function:

logLik(Mlda)

# Let us estimate the log-likelihood for 2,3,4, and 5 cluster solutions: 

lg <- list()

for (i in 2:5) {

Mlda <- LDA(M, k = i, control = list(alpha = 10, delta = .1, seed = 68), method = "Gibbs")

lg[[i]] <- logLik(Mlda) 

    }

plot(2:5, unlist(lg))   

Interpreting Clusters and Dimensions

First, we investigate the relationship between LDA clusters / SVD dimensions and psycho-demographic traits using correlation. 

# Correlate user traits and their SVD scores

# users[,-1] is used to exclude the column with IDs

cor(u_rot, users[,-1], use = "pairwise")

# LDA version

cor(gamma, users[,-1], use = "pairwise")

You can also print the Likes with the highest and lowest varimax-rotated SVD and LDA scores:

# SVD

top <- list()

bottom <-list()

for (i in 1:5) {

f <- order(v_rot[ ,i])

temp <- tail(f, n = 10)

top[[i]]<-colnames(M)[temp]  

temp <- head(f, n = 10)

bottom[[i]]<-colnames(M)[temp]  

}

# LDA

top <- list()

for (i in 1:5) {

f <- order(beta[i,])

temp <- tail(f, n = 10)

top[[i]]<-colnames(M)[temp]  

}

Building Prediction Models

# Split users into 10 groups

folds <- sample(1:10, size = nrow(users), replace = T)

# Take users from group 1 and assign them to the TEST subset

test <- folds == 1

# Extract SVD dimensions from the TRAINING subset

# training set can be accessed using !test

Msvd <- irlba(M[!test, ], nv = 50)

# Rotate Like SVD scores (V)

v_rot <- unclass(varimax(Msvd$v)$loadings)

# Rotate user SVD scores *for the entire sample*

u_rot <- as.data.frame(as.matrix(M %*% v_rot))

# Build linear regression model for openness

# using TRAINING subset

fit_o <- glm(users$ope~., data = u_rot, subset = !test)

# Inspect the regression coefficients

coef(fit_o)

# Do the same for gender

# use family = "binomial" for logistic regression model

fit_g <- glm(users$gender~.,data = u_rot, subset = !test, family = "binomial")

Next, examine the accuracy of the models:

# Compute the predictions for the TEST subset

pred_o <- predict(fit_o, u_rot[test, ])

pred_g <- predict(fit_g, u_rot[test, ], type = "response")

# Correlate predicted and actual values for the TEST subset

r <- cor(users$ope[test], pred_o)

r

# Compute Area Under the Curve for gender

# remember to install ROCR library first

library(ROCR)

temp <- prediction(pred_g, users$gender[test])

auc <- performance(temp,"auc")@y.values

auc

You can also examine the accuracy across different values of k. 

# Choose which k are to be included in the analysis

ks<-c(2:10,15,20,30,40,50)

# Preset an empty list to hold the results

rs <- list()

# Run the code below for each k in ks

for (k in ks){

# Varimax rotate Like SVD dimensions 1 to k

v_rot <- unclass(varimax(Msvd$v[, 1:k])$loadings)

# This code is exactly like the one discussed earlier

u_rot <- as.data.frame(as.matrix(M %*% v_rot))

fit_o <- glm(users$ope~., data = u_rot, subset = !test)

pred_o <- predict(fit_o, u_rot[test, ])

# Save the resulting correlation coefficient as the

# element of R called k

rs[[as.character(k)]] <- cor(users$ope[test], pred_o)

}

# Check the results

rs

Also, let's compute prediction accuracy (for openness) across all 10 cross-validation folds:

pred_o <- rep(NA, n = nrow(users))

for (i in 1:10){

  test <- folds == i

  Msvd <- irlba(M[!test, ], nv = 50)

  v_rot <- unclass(varimax(Msvd$v)$loadings)

  u_rot <- as.data.frame(as.matrix(M %*% v_rot))

  fit_o <- glm(users$ope~., data = u_rot, subset = !test)

  pred_o[test] <- predict(fit_o, u_rot[test, ])

}

r <- cor(users$ope, pred_o)

Finally, the following code can be used to estimate user LDA cluster scores in a cross-validated way:

Mlda <- LDA(M[!test, ], control=list(alpha=1, delta=.1, seed=68), k=50, method="Gibbs")

temp<-posterior(Mlda, M)

gamma<-as.data.frame(temp$topics)

Cross-Validated Predictions of All Traits

Below, you will find a complete code that can be used to build prediction models for all of the traits.

# Load libraries

require(Matrix)

library(ROCR)

library(topicmodels)

library(irlba)

# Load files

users<-read.csv("users.csv")

likes<-read.csv("likes.csv")

ul<-read.csv("users-likes.csv")

# Construct the matrix

ul$user_row<-match(ul$userid,users$userid)

ul$like_row<-match(ul$likeid,likes$likeid)

M<-sparseMatrix(i=ul$user_row,j=ul$like_row,x=1)

rownames(M)<-users$userid

colnames(M)<-likes$name

rm(ul,likes)

# Matrix trimming

while (T){

  i<-sum(dim(M))

  M<-M[rowSums(M)>=50, colSums(M)>=150]

  if (sum(dim(M))==i) break

}

users <- users[match(rownames(M),users$userid), ]

# Start predictions

set.seed(seed=68)

n_folds<-10                # set number of folds

k<-50                      # set k

vars<-colnames(users)[-1]  # choose variables to predict

folds <- sample(1:n_folds, size = nrow(users), replace = T)

results<-list()

for (fold in 1:n_folds){ 

  print(paste("Cross-validated predictions, fold:", fold))

  test <- folds == fold

  # If you want to use SVD:

  Msvd <- irlba(M[!test, ], nv = k)

  v_rot <- unclass(varimax(Msvd$v)$loadings)

  predictors <- as.data.frame(as.matrix(M %*% v_rot))

  # If you want to use LDA, comment out the SVD lines above, and uncomment two lines below

  # Mlda <- LDA(M[!test, ], control = list(alpha = 1, delta = .1, seed=68), k = k, method = "Gibbs")

  # predictors <- as.data.frame(posterior(Mlda,M, control = list(alpha = 1, delta = .1))$topics)

  for (var in vars){

    results[[var]]<-rep(NA, n = nrow(users))

    # check if the variable is dichotomous

    if (length(unique(na.omit(users[,var]))) ==2) {    

      fit <- glm(users[,var]~., data = predictors, subset = !test, family = "binomial")

      results[[var]][test] <- predict(fit, predictors[test, ], type = "response")

    } else {

      fit<-glm(users[,var]~., data = predictors, subset = !test)

      results[[var]][test] <- predict(fit, predictors[test, ])

    }

    print(paste(" Variable", var, "done."))

  }

}

compute_accuracy <- function(ground_truth, predicted){

  if (length(unique(na.omit(ground_truth))) ==2) {

    f<-which(!is.na(ground_truth))

    temp <- prediction(predicted[f], ground_truth[f])

    return(performance(temp,"auc")@y.values)

  } else {return(cor(ground_truth, predicted,use = "pairwise"))}

}

accuracies<-list()

for (var in vars) accuracies[[var]]<-compute_accuracy(users[,var], results[[var]])