Module 2:Preference Modeling with Notes

library("MPsychoR")
data("bandpref")
bandpref

library("BradleyTerry2")

package 㤼㸱BradleyTerry2㤼㸲 was built under R version 4.0.5

bandsBT <- BTm(cbind(Win1, Win2),Band1, Band2, data = bandpref)
bandsAbil <- BTabilities(bandsBT)
round(sort(bandsAbil[,1]), 3)

Scorpions   Emperor    Slayer     Death      Rush 
   -0.312    -0.024     0.000     0.199     0.380 

alphas <- exp(bandsAbil[,1])
sort(alphas)

Scorpions   Emperor    Slayer     Death      Rush 
0.7321361 0.9759494 1.0000000 1.2200734 1.4621035 

alphas["Rush"]/alphas["Death"]

    Rush 
1.198373 

library("psychotree")
data("Topmodel2007")
topmtree <- bttree(preference ~ age + gender + q1 + q2 + q3,
data = Topmodel2007)

plot(topmtree)

Here you can see where the preferences break among questions. The first break is with age 52. Those older than age 52 found Barbra to be the most attractive, and Anni to be the least (this is a terrible example - however, it was the data that was used.) The second split is with q2 - which asks if they watch the show often. There is a final split among males and females. I find this technique and diagram to be super helpful.

pcmat <- as.matrix(Topmodel2007[[1]])
pcvec <- as.vector(pcmat)
pcvec[pcvec == -1] <- 0

library("stringr")
modnames <- t(str_split(colnames(pcmat), ":", simplify = TRUE))

preds <- Topmodel2007[2:6]
ind <- sapply(preds, is.factor)
preds[ind] <- sapply(preds[ind],
function(f) c(0:(length(levels(f))-1))[f])
preds <- scale(preds)
rownames(preds) <- paste0("P", 1:nrow(preds))
head(preds)

       gender        age         q1         q2         q3
P1 -0.9973924  2.1396096  1.1675757  0.4709234  0.8796174
P2 -0.9973924 -0.9476026 -0.8520147 -2.1124278 -1.1309367
P3  0.9973924 -0.9476026 -0.8520147 -2.1124278 -1.1309367
P4  0.9973924 -1.0162074 -0.8520147  0.4709234 -1.1309367
P5 -0.9973924 -0.9476026 -0.8520147 -2.1124278 -1.1309367
P6  0.9973924 -1.0162074 -0.8520147  0.4709234 -1.1309367

library("BTLLasso")
sid <- rep(rownames(preds), ncol(pcmat))
mfirst <- rep(modnames[1,], each = nrow(preds))
msecond <- rep(modnames[2,], each = nrow(preds))
BTresp <- response.BTLLasso(response = pcvec, subject = sid,
first.object = mfirst, second.object = msecond)

set.seed(123)
lambda <- exp(seq(log(10), log(1), length = 20))-1
modLasso <- cv.BTLLasso(Y = BTresp, X = preds, lambda = lambda,
folds = 5, trace = FALSE)

Full model 
Cross-Validation... 

modLasso

Output of BTLLasso estimation: 
--- 
Setting:
 192 subjects
 6 objects
 2 response categories
 5 subject-specific covariate(s)
 0 subject-object-specific covariate(s) with object-specific effects
 0 (subject-)object-specific covariate(s) with global effects
 No order effect
 20 different tuning parameters 

 Cross-validation criterion: RPS 
--- 
Parameter estimates after 5 - fold cross-validation 

Intercepts: 
   Anja    Anni Barbara   Fiona    Hana   Mandy 
 -0.372  -0.097   0.351   0.185   0.402  -0.468 

Object-specific effects for subject-specific covariate(s): 
         Anja   Anni Barbara  Fiona   Hana  Mandy
gender -0.119  0.207   0.207 -0.057 -0.119 -0.119
age     0.166 -0.329   0.166 -0.142 -0.026  0.166
q1      0.177  0.131   0.016  0.016 -0.222 -0.119
q2      0.239 -0.105   0.086  0.171 -0.286 -0.105
q3     -0.166  0.031  -0.166  0.031  0.136  0.136

--- 

Optimal lambda: 0.6237767 

Log likelihood: -869.9134 

plot(modLasso)

optind <- which.min(modLasso$criterion)
lambda2 <- lambda[(optind-3):(optind+3)]
set.seed(111)
bootLasso <- boot.BTLLasso(modLasso, lambda = lambda2,
cores = 4, B = 50, trace = FALSE, trace.cv = FALSE)

plot(bootLasso, plots_per_page = 6, ask_new = FALSE)

These are confidence interval, which are calcuated using a bootstrap method.