Minecraft_Markdown_FINAL.Rmd

---
title: "Minecraft"
author: "Heinrich Peters"
date: "11/8/2019"
output: html_document
---

```{r setup, include=FALSE}
# Set knitting options
knitr::opts_chunk$set(echo = TRUE)

# MAC
 knitr::opts_knit$set(root.dir = '/Users/hp2500/Google Drive/STUDY/Columbia/Research/Minecraft/Data')

#PC
#knitr::opts_knit$set(root.dir = 'C:\\Users\\Heinrich\\Google Drive\\STUDY\\Columbia\\Research\\Minecraft\\Data')

# Set graphic parameters
par(ask=T)

```


# Import Libraries
```{r, results = "hide", message=FALSE}

library(plyr)
library(reshape2)
library(ltm)
library(class)
library(pscl)
library(paran)
library(eRm)
library(foreign)
library(xtable)
library(stargazer)
library(plyr)
library(psych)
library(rlang)
library(ggplot2)
library(mirt)
library(lavaan)
library(semPlot)
library(colorspace)
library(missForest)
library(effsize)
library(glmnet)
library(BaylorEdPsych)
library(mvnmle)
library(MissMech)
library(normtest)
library(cocor)
library(CTT)
library(MVN)
library(tidyverse)
library(Hmisc) 
library(mlr3)
library(mlr3viz)
library(mlr3learners)
library(mlr3tuning)
library(mlr3misc)



```


# Define functions
```{r}

# Define functions to aggregate data
aggdata <- function(dataset){
  ddply(dataset, c("id","item"), summarise, 
        n = length(id),
        time = max(as.numeric(time_taken)),
        distance = max(as.numeric(distance_travelled)),
        answer = (answer[n]) 
  )
}

aggdata_dem <- function(dataset2){
  ddply(dataset2, c("id"), summarise,
        age = max(Age),
        gender = max(Gender),
        gamesGen = max(gamesGen),
        gamingSkills = max(gamingSkills),
        mcExperience = max(mcExperience),
        gameFun = max(gameFun),
        takeAgain = max(takeAgain),
        Recommend = max(Recommend),
        Distracted = max(Distracted),
        Boring = max(Boring)
  )
}


# Define function to convert long to wide format
long_to_wide <- function(data, x){
  wide <- dcast(data, id~item, value.var = "answer")
  #wide[is.na(wide)]<-"Timeout"
  
  wide[wide=="Timeout"]<-0
  wide[wide=="False"]<-0
  wide[wide=="True"]<-1
  
  names(wide) <- paste0(c("",rep(x,12)),names(wide))
  wide <- as.data.frame(sapply(wide, as.numeric))
  return(wide)
}

# Define function for correlation tables
corstarsl <- function(x){ 
  x <- as.matrix(x) 
  R <- rcorr(x, type ="spearman")$r 
  p <- rcorr(x, type ="spearman")$P 
  
  # define notions for significance levels; spacing is important.
  mystars <- ifelse(p < .001, "***", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " ")))
  
  # trunctuate the matrix that holds the correlations to two decimal
  R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1] 
  
  # build a new matrix that includes the correlations with their apropriate stars 
  Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x)) 
  diag(Rnew) <- paste(diag(R), " ", sep="") 
  rownames(Rnew) <- colnames(x) 
  colnames(Rnew) <- paste(colnames(x), "", sep="") 
  
  # remove upper triangle
  Rnew <- as.matrix(Rnew)
  Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
  Rnew <- as.data.frame(Rnew) 
  
  # remove last column and return the matrix (which is now a data frame)
  Rnew <- cbind(Rnew[1:length(Rnew)-1])
  return(Rnew) 
}

```


# CLEAN PC DATA
```{r}
# Read pattern completion data
# Get the files names
files_pc = list.files(pattern="*Pattern")

# First apply read.csv, then rbind
pc = do.call(rbind, lapply(files_pc, function(x) read.csv(x, sep=",", stringsAsFactors = F)))

# aggregate data to item level
pc_agg <- aggdata(pc)
head(pc_agg)

# Convert long to wide format
pc_wide <- long_to_wide(pc_agg, "pc")
head(pc_wide)

```

# CLEAN MR DATA
```{r}
# Read mental rotation data
# Get the files names mr
files_mr = list.files(pattern="*Mental")

# First apply read.csv, then rbind
mr = do.call(rbind, lapply(files_mr, function(x) read.csv(x, sep=",", stringsAsFactors = F)))

# aggregate data to item level
mr_agg <- aggdata(mr)
head(mr_agg)

# Identify cashed items
mr_agg[mr_agg$n < 20, ]$answer <- NA

# Convert long to wide format
mr_wide <- long_to_wide(mr_agg, "mr")
head(mr_wide)
```

# CLEAN SC DATA
```{r, warning=FALSE}
# Read structure recreation data
# Get the files names sc
files_sc = list.files(pattern="*Rec")

# First apply read.csv, then rbind
sc = do.call(rbind, lapply(files_sc, function(x) read.csv(x, sep=",", stringsAsFactors = F)))

# aggregate data to item level
sc_agg <- aggdata(sc)
head(sc_agg)

# Identify cashed items
sc_agg[sc_agg$n < 20, ]$answer <- NA

# Convert long to wide format
sc_wide <- long_to_wide(sc_agg, "sc")
head(sc_wide)
```

# CREATE DATASETS
```{r}

# Merge to ensure same order of subjects in all data frames
all_wide <- merge(merge(pc_wide, mr_wide, by = 'id'), sc_wide, by = 'id')
all_wide <- all_wide %>% map(as_factor) %>% as.data.frame()

pc_wide <- all_wide[,1:13]
mr_wide <- all_wide[,c(1,14:25)]
sc_wide <- all_wide[,c(1,26:37)]

# Percentage missing
sum(is.na(all_wide))/prod(dim(all_wide))
sum(is.na(pc_wide))/prod(dim(pc_wide))
sum(is.na(mr_wide))/prod(dim(mr_wide))
sum(is.na(sc_wide))/prod(dim(sc_wide))

# MCAR Test
LittleMCAR(pc_wide)$p
LittleMCAR(mr_wide)$p
LittleMCAR(sc_wide)$p

# impute missing data
set.seed(1)
all_imp <- missForest(all_wide[,-1],ntree = 1000)$ximp
all_wide <- cbind(all_wide[,1], all_imp)
all_wide <- as.data.frame(lapply(all_wide, function(x) as.numeric(as.character(x))))
names(all_wide)[1] <- "id"
head(all_wide)

# split data frame into subtest level data frames
pc_wide <- all_wide[,1:13]
mr_wide <- all_wide[,c(1,14:25)]
sc_wide <- all_wide[,c(1,26:37)]

names(pc_wide)[1] <- "id"
names(mr_wide)[1] <- "id"
names(sc_wide)[1] <- "id"

# Create vectors for differential item functioning
dem <- aggdata_dem(pc)
gender <- merge(all_wide,dem)[,c(1,39)]
mcexp <- merge(all_wide,dem)[,c(1,42)]

```


# ITEM ANALYSIS PC (eRm)

```{r}
#### Exclude Items ####
# training items 
pc_wide <- subset(pc_wide, select= -c(pc1))
pc_wide <- subset(pc_wide, select= -c(pc2)) 

# low variance
pc_wide <- subset(pc_wide, select= -c(pc4)) 

# low item fit
pc_wide <- subset(pc_wide, select= -c(pc12))

#### Descriptive stats ####
descript_pc <- descript(pc_wide[,-1])
descript_pc
summary(descript_pc$ExBisCorr)
sd(descript_pc$ExBisCorr)

#### Fit Rasch model ####
# Item parameters
res_pc <- RM(pc_wide[,-1])
res_pc
-coef(res_pc)
summary(res_pc)

# Person parameters
pres_pc <- person.parameter(res_pc)
pc_pers <- cbind(pc_wide$id, coef(pres_pc))

# Stepwise item elimination
stepwiseIt(res_pc, criterion = list("itemfit"), alpha = 0.05,
           verbose = TRUE, maxstep = NA)

mean(pc_pers[,-1])
sd(pc_pers[,-1])

#### Model tests ####
# Item Fit // Person Fit
itemfit_pc <- eRm::itemfit(pres_pc)
itemfit_pc
personfit_pc <- eRm::personfit(pres_pc)

# Wald Test
Waldtest(res_pc)

# Goodness of Fit
gofIRT(pres_pc)

# Andersen Likelihood Ratio Test 
lrres_pc<- LRtest(res_pc, splitcr = "mean")
lrres_pc
lrres_pc<- LRtest(res_pc, splitcr = "median")
lrres_pc


# Follow up wald test
Waldtest(res_pc, splitcr = gender[,2])

# Informaion Criteria
IC(pres_pc)

#### Visualization / plots ####
plotjointICC(res_pc, main="Pattern Completion", lwd=1, xlim= c(-5,5))
plotINFO(res_pc, legpos=c(10,10), type = "both")


#### Test for unidimensionality ####
paran(pc_wide[,-1])  # is that ok?
pc_unidim <- rasch(pc_wide[,-1])
unidimTest(pc_unidim)

#### Prepare for LaTex export ####
p_value_pc <- as.vector(round(sapply(itemfit_pc$i.fit, function(x) 1-pchisq(x,116)),3))
pc_ia_tab <- data.frame(-res_pc$betapar,res_pc$se.beta,descript_pc$ExBisCorr,itemfit_pc$i.fit, itemfit_pc$i.df-1, p_value_pc,  row.names = NULL)
pc_ia_tab <- round(pc_ia_tab, 3)
pc_ia_tab <- data.frame(names(pc_wide[,-1]), pc_ia_tab)
names(pc_ia_tab) <- c("item", "beta", "se", "cor", "chisq", "df", "p-value")
pc_ia_tab <- pc_ia_tab[order(pc_ia_tab[,2]),]
xtable(data.frame(pc_ia_tab, row.names =NULL), type = "latex")

```

# ITEM ANALYSIS MR (eRm)
```{r}
#### Exclude items ####
# training items
mr_wide <- subset(mr_wide, select=-c(mr1))  
mr_wide <- subset(mr_wide, select=-c(mr2))

# low item fit
mr_wide <- subset(mr_wide, select=-c(mr9))  
mr_wide <- subset(mr_wide, select=-c(mr12)) 
mr_wide <- subset(mr_wide, select=-c(mr11)) 

#### Descriptive stats ####
descript_mr <- descript(mr_wide[,-1])
descript_mr
mean(descript_mr$ExBisCorr)
sd(descript_mr$ExBisCorr)
summary(descript_mr$ExBisCorr)

#### Fit Rasch model ####
# Item parameters
res_mr <- RM(mr_wide[,-1])
res_mr
-coef(res_mr)
summary(res_mr)

# Person parameters
pres_mr <- person.parameter(res_mr)
mr_pers <- cbind(mr_wide$id, coef(pres_mr))

# Stepwise item elimination
stepwiseIt(res_mr, criterion = list("itemfit"), alpha = 0.05,
           verbose = TRUE, maxstep = NA)

mean(mr_pers[,-1])
sd(mr_pers[,-1])

#### Model fit ####

# Item Fit // Person Fit
itemfit_mr <- eRm::itemfit(pres_mr)
personfit_mr <- eRm::personfit(pres_mr)

# Wald Test
Waldtest(res_mr)

# Goodness of Fit
gofIRT(pres_mr)

# Andersen Likelihood Ratio Test 
lrres_mr<- LRtest(res_mr, splitcr = "mean")
lrres_mr
lrres_mr<- LRtest(res_mr, splitcr = "median")
lrres_mr

# Informaion Criteria
IC(pres_mr)


#### Plot / visualization ####
#plotICC(res_mr)
plotjointICC(res_mr, legend = T, main="Mental Rotation", lwd=2)
plotINFO(res_mr, legpos=c(10,10), type = "both")


#### Test for unidimensionality ####
paran(mr_wide[,-1])
fit_mr <- rasch(mr_wide[,-1])
unidimTest(fit_mr)

#### Prepare for LaTex export ####
p_value_mr <- as.vector(round(sapply(itemfit_mr$i.fit, function(x) 1-pchisq(x,itemfit_mr$i.df[1]-1)),3))
mr_ia_tab <- data.frame(-res_mr$betapar,res_mr$se.beta,descript_mr$ExBisCorr,itemfit_mr$i.fit, itemfit_mr$i.df-1, p_value_mr,  row.names = NULL)
mr_ia_tab <- round(mr_ia_tab, 3)
mr_ia_tab <- data.frame(names(mr_wide[,-1]), mr_ia_tab)
names(mr_ia_tab) <- c("item", "beta", "se", "cor", "chisq", "df", "p-value")
mr_ia_tab <- mr_ia_tab[order(mr_ia_tab[,2]),]
xtable(mr_ia_tab)

```

# ITEM ANALYSIS SC (eRm)
```{r}
#### Exclude Items ####
# Training items first
sc_wide <- subset(sc_wide, select= -c(sc1))
sc_wide <- subset(sc_wide, select= -c(sc2))

# Low variance
sc_wide <- subset(sc_wide, select= -c(sc3))
sc_wide <- subset(sc_wide, select= -c(sc5))

# Low item fit
sc_wide <- subset(sc_wide, select= -c(sc11))
sc_wide <- subset(sc_wide, select= -c(sc4))

#### Descriptive stats ####
descript_sc <- descript(sc_wide[,-1])
descript_sc
mean(descript_sc$ExBisCorr)
sd(descript_sc$ExBisCorr)
summary(descript_sc$ExBisCorr)

#### Fit Rasch model ####
# Item parameters
res_sc <- RM(sc_wide[,-1])
res_sc
-coef(res_sc)
summary(res_sc)

# Person parameters
pres_sc <- person.parameter(res_sc)
sc_pers <- data.frame(sc_wide$id, coef(pres_sc))

# Stepwise item elimination 
stepwiseIt(res_sc, criterion = list("itemfit"), alpha = 0.05,
           verbose = TRUE, maxstep = NA)

mean(sc_pers[,-1])
sd(sc_pers[,-1])


#### Model fit ####
# Item Fit // Person Fit
itemfit_sc <- eRm::itemfit(pres_sc)
itemfit_sc
personfit_sc <- eRm::personfit(pres_sc)

# Goodness of Fit
gofIRT(pres_sc)

# Andersen Likelihood Ratio Test 
lrres_sc<- LRtest(res_sc, splitcr = "mean")
lrres_sc
lrres_sc<- LRtest(res_sc, splitcr = "median")
lrres_sc



# Wald Test
Waldtest(res_sc, splitcr="mean")
Waldtest(res_sc, splitcr="median")

# Informaion Criteria
IC(pres_sc)

#### Plot / visualization ####
#plotICC(res_sc)
plotjointICC(res_sc, legend =T, xlim = c(-6,6), main = "Structure Recreation", col=c("black", "red", "blue", "green", "magenta", "purple", "orange"), lwd = 2)
plotINFO(res_sc, legpos=c(10,10), type = "both")

#### Test for unidimensionality ####
paran(sc_wide[,-1])
fit_sc <- rasch(sc_wide[,-1])
fit_sc
unidimTest(fit_sc)
dev.off()

#### Prepare for LaTex export ####
p_value_sc <- as.vector(round(sapply(itemfit_sc$i.fit, function(x) 1-pchisq(x,itemfit_sc$i.df[1]-1)),3))
sc_ia_tab <- data.frame(-res_sc$betapar,res_sc$se.beta,descript_sc$ExBisCorr,itemfit_sc$i.fit, as.integer(itemfit_sc$i.df-1), p_value_sc, row.names = NULL)
sc_ia_tab <- round(sc_ia_tab, 3)
sc_ia_tab <- data.frame(names(sc_wide[,-1]), sc_ia_tab)
names(sc_ia_tab) <- c("item", "beta", "se", "cor", "chisq", "df", "p-value")
sc_ia_tab <- sc_ia_tab[order(sc_ia_tab[,2]),]
xtable(sc_ia_tab)
```


# FACTOIAL VALIDITY
```{r}
#### Confirmatory factor analysis ####
all_wide_2 <- merge(merge(pc_wide, mr_wide), sc_wide)
head(all_wide_2)
dim(all_wide_2)

# 3 factor solution oblique
cfa_3f <- 'pc =~ pc3 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11
mr =~ mr3 + mr4 + mr5 + mr6 + mr7 + mr8 + mr10
sc =~ sc6 + sc7 +sc8 + sc9+ sc10 + sc12'
fit_cfa_3f <- cfa(cfa_3f, data = all_wide_2[,-1], std.lv=T, orthogonal = F)
summary(fit_cfa_3f, standardized = T, fit.measures = T, rsq = T)
semPaths(fit_cfa_3f, "std",cardinal = F, edge.label.cex = 0.6, curvePivot =T, layout = "tree2", sizeMan = 3, sizeMan2 = 2, edge.color = "black", edge.label.position = 0.8)

# 3 factor solution orthogonal
fit_cfa.ort <- cfa(cfa_3f, data = all_wide_2[,-1], std.lv=T, orthogonal = T)
summary(fit_cfa.ort, standardized = T, fit.measures = T, rsq = T)
anova(fit_cfa_3f, fit_cfa.ort) 

# compare to model with just one latent factor
cfa.model.one <- 'score =~ pc3 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11 + mr3 + mr4 + mr5 + mr6 + mr7 + mr8 + mr10 + sc6 + sc7 +sc8 + sc9 + sc10 + sc12'
fit_cfa.one <- cfa(cfa.model.one, data = all_wide_2[,-1])
summary(fit_cfa.one, standardized = T, fit.measures = T, rsq = T)
anova(fit_cfa_3f, fit_cfa.one)


# 2 factor solution oblique MR excluded
cfa_2f <- 'pc =~ pc3 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11
sc =~ sc6 + sc7 +sc8 + sc9+ sc10 + sc12'
fit_cfa <- cfa(cfa_2f, data = all_wide_2[,-1], std.lv=T, orthogonal = F)
summary(fit_cfa, standardized = T, fit.measures = T, rsq = T)
semPaths(fit_cfa, "std",cardinal = F, edge.label.cex = 0.6, curvePivot =T, layout = "tree2", sizeMan = 3, sizeMan2 = 2, edge.color = "black", edge.label.position = 0.8)

# 2 factor solution orthogonal MR excluded
cfa_2f_orth <- 'pc =~ pc3 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11
sc =~ sc6 + sc7 +sc8 + sc9+ sc10 + sc12'
fit_cfa_orth <- cfa(cfa_2f_orth, data = all_wide_2[,-1], std.lv=T, orthogonal = T)
summary(fit_cfa_orth, standardized = T, fit.measures = T, rsq = T)
semPaths(fit_cfa_orth, "std",cardinal = F, edge.label.cex = 0.6, curvePivot =T, layout = "tree2", sizeMan = 3, sizeMan2 = 2, edge.color = "black", edge.label.position = 0.8)
anova(fit_cfa, fit_cfa_orth)

# compare to model with just one latent factor
cfa.model.one <- 'score =~ pc3 + pc5 + pc6 + pc7 + pc8 + pc9 + pc10 + pc11 + sc6 + sc7 +sc8 + sc9 + sc10 + sc12'
fit_cfa.one <- cfa(cfa.model.one, data = all_wide_2[,-1])
summary(fit_cfa.one, standardized = T, fit.measures = T, rsq = T)
anova(fit_cfa, fit_cfa.one)
```

# CONVERGENT VALIDITY
```{r, warning=FALSE}
#### Prepare data ####
#Read SPM and VKMR data
# dat_pen_paper <- read.csv("C:\\Users\\Heinrich\\Google Drive\\DATA-20180216T104535Z-001\\dat_merged")
dat_pen_paper <- read.csv("dat_merged")
dat_pen_paper$X <- NULL
names(dat_pen_paper)[1] <- "id"

# some descriptive stats of paper based tests
mean(dat_pen_paper$VKMR_scores)
sd(dat_pen_paper$VKMR_scores)
mean(dat_pen_paper$SPM_scores)
sd(dat_pen_paper$SPM_scores)

# Prepare MARS data for merging
pc_pers <- data.frame(pc_wide$id, coef(pres_pc))
mr_pers <- data.frame(mr_wide$id, coef(pres_mr))
sc_pers <- data.frame(sc_wide$id, coef(pres_sc))

names(pc_pers)[1] <- "id"
names(mr_pers)[1] <- "id"
names(sc_pers)[1] <- "id"

# Merge MARS data with SPM and VKMR data
df_pers<- join_all(list(pc_pers, mr_pers, sc_pers), by = 'id')
df_all<- join_all(list(df_pers,dat_pen_paper), by = 'id', type = "inner")
names(df_all) <- c("id", "pc", "mr", "sc", "VKMR", "SPM")
head(df_all)
dim(df_all)

#### Check assumptions ####
# multivariate and univariate normality
mvn(df_all[,-1], mvnTest = "hz")

#### Correlation analysis ####
# Correlation matrices
round(cor(df_all[,-1], use = "pairwise.complete"),2)
round(cor(df_all[,-1], use = "pairwise.complete", method = "spearman"),2)
corstarsl(df_all[,-1]) # significance levels for spearman correlations

# Correlation matrix with attenuation correction
rel <- c(0.74, 0.61, 0.76, 0.90, 0.76)
round(correct.cor(cor(df_all[,-1]), rel),2)
round(correct.cor(cor(df_all[,-1], method = "spearman"), rel),2)

# Export for LaTex
xtable(round(correct.cor(cor(df_all[,-1]), rel),2))
xtable(round(correct.cor(cor(df_all[,-1], method = "spearman"), rel),2))

# Significance tests
cor.test(df_all$pc,df_all$SPM, method = "spearman", exact=F, alternative = "greater")
cor.test(df_all$pc,df_all$VKMR, method = "spearman", exact=F, alternative = "greater")
cor.test(df_all$mr,df_all$VKMR, method = "spearman", exact=F, alternative = "greater")
cor.test(df_all$mr,df_all$SPM, method = "spearman", exact=F, alternative = "greater")
cor.test(df_all$sc,df_all$VKMR, method = "spearman", exact=F, alternative = "greater")
cor.test(df_all$sc,df_all$SPM, method = "spearman", exact=F, alternative = "greater")

#### Structural equation model #### 
SEM_2f <- 'MARS =~ pc + sc
          PB =~ VKMR + SPM
          MARS ~ PB'

fit_SEM_2f <- sem(SEM_2f, data = df_all[,-1], std.lv=F, orthogonal = F, estimator = "MLM")  # MLM because not multivariate nomal
summary(fit_SEM_2f, standardized = T, fit.measures = T, rsq = T)
semPaths(fit_SEM_2f, "std",cardinal = F, edge.label.cex = 0.6, curvePivot =T, layout = "tree", sizeMan = 3, sizeMan2 = 2, edge.color = "black", edge.label.position = 0.8)

```




# Discriminant Validity
```{r}
# compare criterion correlations pc
r.jk <- cor(df_all$pc,df_all$SPM, method = "spearman")
r.jk
r.jh <- cor(df_all$pc,df_all$VKMR, method = "spearman")
r.jh
r.kh <- cor(df_all$SPM,df_all$VKMR, method = "spearman")
r.kh
cocor.dep.groups.overlap(r.jk, r.jh, r.kh, n = 120, alternative = "greater")

# compare criterion correlations mr
r.jk <- cor(df_all$mr,df_all$VKMR, method = "spearman")
r.jk
r.jh <- cor(df_all$mr,df_all$SPM, method = "spearman")
r.jh
r.kh <- cor(df_all$SPM,df_all$VKMR, method = "spearman")
r.kh
cocor.dep.groups.overlap(r.jk, r.jh, r.kh, n = 120, alternative = "greater")

# compare criterion correlations sc
r.jk <- cor(df_all$sc,df_all$VKMR, method = "spearman")
r.jk
r.jh <- cor(df_all$sc,df_all$SPM, method = "spearman")
r.jh
r.kh <- cor(df_all$SPM,df_all$VKMR, method = "spearman")
r.kh
cocor.dep.groups.overlap(r.jk, r.jh, r.kh, n = 120, alternative = "greater")
```

# Demographics, experience and enjoyment
```{r}
# ceate new dataset with experience and demographic variables
teq <- aggdata_dem(pc)
head(teq)
summary(teq)
dim(teq)

summary(teq$age)
sd(teq$age)
table(teq$gender)
table(teq$mcExperience)
prop.table(table(teq$mcExperience))

summary(teq$gamesGen)
summary(teq$gamingSkills)

df_all_teq <- merge(df_all, teq)
head(df_all_teq)
dim(df_all_teq)

summary(df_all_teq$age)
sd(df_all_teq$age)
table(df_all_teq$gender)
table(df_all_teq$mcExperience)
prop.table(table(df_all_teq$mcExperience))

``` 

# Item analysis test enjoyment questionnaire
```{r}
# shrink dataset to relevant variables
teq_2 <- subset(teq, select = -c(age,gender, gamesGen, gamingSkills,mcExperience))
head(teq_2)

# recode negatively coded variables
teq_2$Distracted <- 6-teq_2$Distracted
teq_2$Boring <- 6-teq_2$Boring
head(teq_2)

teq_agg <- data.frame(teq_2$id,rowSums(teq_2[,-1]))
names(teq_agg) <- c("id","liking")
head(teq_agg)

hist(teq_agg[,2], main = "Test Enjoyment Scale", xlab = "Scores")

df_all_teq2 <- merge(df_all_teq, teq_agg) 
head(df_all_teq2)#

mean(df_all_teq2$liking)
sd(df_all_teq2$liking)


```


# Test for effects of gender
```{r}
bonf = 8

# pc~gender
t_pc_gender = t.test(df_all_teq$pc~df_all_teq$gender)
t_pc_gender
wilcox.test(df_all_teq$pc~df_all_teq$gender)
cohen.d(df_all_teq$pc~factor(df_all_teq$gender))
t_pc_gender$p.value*bonf

# mr~gender
t_mr_gender = t.test(df_all_teq$mr~df_all_teq$gender, alternative = "less")
t_mr_gender
wilcox.test(df_all_teq$mr~df_all_teq$gender, alternative = "less")
cohen.d(df_all_teq$mr~factor(df_all_teq$gender))
t_mr_gender$p.value*bonf

# sc~gender
t_sc_gender = t.test(df_all_teq$sc~df_all_teq$gender, alternative = "less")
t_sc_gender 
wilcox.test(df_all_teq$sc~df_all_teq$gender, alternative = "less")
cohen.d(df_all_teq$sc~factor(df_all_teq$gender))
t_sc_gender$p.value*bonf

# spm~gender
t_spm_gender = t.test(df_all_teq$SPM~df_all_teq$gender)
t_spm_gender
wilcox.test(df_all_teq$SPM~df_all_teq$gender)
cohen.d(df_all_teq$SPM~factor(df_all_teq$gender))
t_spm_gender$p.value*bonf

# vkmr~gender
t_vkmr_gender = t.test(df_all_teq$VKMR~df_all_teq$gender, alternative = "less")
t_vkmr_gender
wilcox.test(df_all_teq$VKMR~df_all_teq$gender, alternative = "less")
cohen.d(df_all_teq$VKMR~factor(df_all_teq$gender))
t_vkmr_gender$p.value*bonf

# enjoyment~gender
t_enj_gender = t.test(data = df_all_teq2, liking~gender)
t_enj_gender
wilcox.test(data = df_all_teq2, liking~gender)
cohen.d(data = df_all_teq2, liking~factor(gender))
t_enj_gender$p.value * bonf

# patg~gender
t_patg_gender = t.test(data = df_all_teq2, gamesGen~gender)
t_patg_gender
wilcox.test(data = df_all_teq2, gamesGen~gender)
cohen.d(data = df_all_teq2, gamesGen~factor(gender))
t_patg_gender$p.value * bonf

# sags~gender
t_sags_gender = t.test(data = df_all_teq2, gamingSkills~gender)
t_sags_gender
wilcox.test(data = df_all_teq2, gamingSkills~gender)
cohen.d(data = df_all_teq2, gamingSkills~factor(gender))
t_sags_gender$p.value *bonf

```

# Test for effects of Minecraft experience
```{r}
bonf = 8

# pc~mcexp
t_pc_mcexp = t.test(df_all_teq$pc~df_all_teq$mcExperience)
t_pc_mcexp
wilcox.test(df_all_teq$pc~df_all_teq$mcExperience)
cohen.d(df_all_teq$pc~factor(df_all_teq$mcExperience))
t_pc_mcexp$p.value*bonf

# mr~mcexp
t_mr_mcexp = t.test(df_all_teq$mr~df_all_teq$mcExperience)
t_mr_mcexp
wilcox.test(df_all_teq$mr~df_all_teq$mcExperience)
cohen.d(df_all_teq$mr~factor(df_all_teq$mcExperience))
t_mr_mcexp$p.value*bonf

# sc~mcexp
t_sc_mcexp = t.test(df_all_teq$sc~df_all_teq$mcExperience, alternative = "less")
t_sc_mcexp
wilcox.test(df_all_teq$sc~df_all_teq$mcExperience, alternative = "less")
cohen.d(df_all_teq$sc~factor(df_all_teq$mcExperience))
t_sc_mcexp$p.value*bonf

# spm~mcexp
t_smp_mcexp = t.test(df_all_teq$SPM~df_all_teq$mcExperience)
t_smp_mcexp
wilcox.test(df_all_teq$SPM~df_all_teq$mcExperience)
cohen.d(df_all_teq$SPM~factor(df_all_teq$mcExperience))
t_smp_mcexp$p.value*bonf

# vkmr~mcexp
t_vkmr_mcexp = t.test(df_all_teq$VKMR~df_all_teq$mcExperience)
t_vkmr_mcexp
wilcox.test(df_all_teq$VKMR~df_all_teq$mcExperience)
cohen.d(df_all_teq$VKMR~factor(df_all_teq$mcExperience))
t_vkmr_mcexp$p.value*bonf

# enj~mcexp
t_enj_mcexp = t.test(data = df_all_teq2, liking~mcExperience)
t_enj_mcexp
wilcox.test(data = df_all_teq2, liking~mcExperience)
cohen.d(data = df_all_teq2, liking~factor(mcExperience))
t_enj_mcexp$p.value*bonf

# patg~mcexp
t_patg_mcexp = t.test(data = df_all_teq2, gamesGen~mcExperience)
t_patg_mcexp
wilcox.test(data = df_all_teq2, gamesGen~mcExperience)
cohen.d(data = df_all_teq2, gamesGen~factor(mcExperience))
t_patg_mcexp$p.value*bonf

# sags~mcexp
t_sags_mcexp = t.test(data = df_all_teq2, gamingSkills~mcExperience)
t_sags_mcexp
wilcox.test(data = df_all_teq2, gamingSkills~mcExperience)
cohen.d(data = df_all_teq2, gamingSkills~factor(mcExperience))
t_sags_mcexp$p.value*bonf

# use one sided t test
```


# Effects of gaming skills and positive attitudes towards games 
```{r}
# Analysis of self assessed gaming skills
cor.test(df_all_teq$pc,df_all_teq$gamingSkills, exact = F, method = "spearman")
cor.test(df_all_teq$mr,df_all_teq$gamingSkills, exact = F, method = "spearman")
cor.test(df_all_teq$sc,df_all_teq$gamingSkills, exact = F, method = "spearman")
cor.test(df_all_teq$SPM,df_all_teq$gamingSkills, exact = F, method = "spearman")
cor.test(df_all_teq$VKMR,df_all_teq$gamingSkills, exact = F, method = "spearman")

# Analysis of positive attitude towards games 
cor.test(df_all_teq$pc,df_all_teq$gamesGen, exact = F, method = "spearman")
cor.test(df_all_teq$mr,df_all_teq$gamesGen, exact = F, method = "spearman")
cor.test(df_all_teq$sc,df_all_teq$gamesGen, exact = F, method = "spearman")
cor.test(df_all_teq$SPM,df_all_teq$gamesGen, exact = F, method = "spearman")
cor.test(df_all_teq$VKMR,df_all_teq$gamesGen, exact = F, method = "spearman")

```

# Group differences after correcting for self assessed gaming skills
```{r}
lm1 <- lm(data = df_all_teq2, sc ~ gamingSkills)


summary(lm1)

t.test(lm1$residuals~ df_all_teq2$gender)
wilcox.test(lm1$residuals~ df_all_teq2$gender)
cohen.d(lm1$residuals~ factor(df_all_teq2$gender))
t.test(lm1$residuals~ df_all_teq2$mcExperience)
wilcox.test(lm1$residuals~ df_all_teq2$mcExperience)
cohen.d(lm1$residuals~ factor(df_all_teq2$mcExperience))
```

# Group differences after correcting for actual minecraft skills
```{r}

# Get times spent on tutorial test item 
files_test = list.files(pattern="*Test")
test = do.call(rbind, lapply(files_test, function(x) read.csv(x, sep=",", stringsAsFactors = F)))
test_agg <- aggdata(test)


df_all_teq_3 = merge(test_agg, df_all_teq2)

lm2 <- lm(data = df_all_teq_3, sc ~ time)
summary(lm2)

cor.test(df_all_teq_3$sc,df_all_teq_3$time)

t.test(lm2$residuals~ df_all_teq_3$gender)
wilcox.test(lm2$residuals~ df_all_teq_3$gender)
cohen.d(lm2$residuals~ factor(df_all_teq_3$gender))
t.test(lm2$residuals~ df_all_teq_3$mcExperience)
wilcox.test(lm2$residuals~ df_all_teq_3$mcExperience)
cohen.d(lm2$residuals~ factor(df_all_teq_3$mcExperience))
```

# Corrected reliabilities 
```{r}
pc_length = 8
mr_length = 7
sc_length = 6

# Number of items needed to reach a reliability of .9
a <- as.numeric(spearman.brown(0.73, 0.9, "r"))
a*pc_length
b <- as.numeric(spearman.brown(0.62, 0.9, "r"))
b*mr_length
c <- as.numeric(spearman.brown(0.76, 0.9, "r"))
c*sc_length

# Reliabilities if scale length was comparable to Quiroga and Foroughi
as.numeric(spearman.brown(0.73, 15/pc_length, "n"))
as.numeric(spearman.brown(0.62, 15/mr_length, "n"))
as.numeric(spearman.brown(0.76, 15/sc_length, "n"))

as.numeric(spearman.brown(0.73, 41/pc_length, "n"))
as.numeric(spearman.brown(0.62, 41/mr_length, "n"))
as.numeric(spearman.brown(0.76, 41/sc_length, "n"))

```

# Additional Plots

```{r}
gg <- ggplot(df_all) + 
  geom_point(mapping = aes(x = VKMR, y = pc)) +
  ggtitle("Relationship between VKMR and PC")

gg

```

```{r}
gg <- ggplot(df_all) + geom_point(mapping = aes(x = VKMR, y = mr)) +
  ggtitle("Relationship between VKMR and MR")

gg

```

```{r}
gg <- ggplot(df_all) + geom_point(mapping = aes(x = VKMR, y = sc)) +
    ggtitle("Relationship between VKMR and SC")

gg

```

```{r}
gg <- ggplot(df_all) + geom_point(mapping = aes(x = SPM, y = pc)) +
    ggtitle("Relationship between VKMR and PC")

gg

```

```{r}
gg <- ggplot(df_all) + geom_point(mapping = aes(x = SPM, y = mr)) +
    ggtitle("Relationship between VKMR and MR")

gg

```

```{r}
gg <- ggplot(df_all) + geom_point(mapping = aes(x = SPM, y = sc)) +
    ggtitle("Relationship between VKMR and SC")

gg

```

# Analyze Process Data 

## Failure to engage 


## Predict Score 
```{r}

vars = c('id', 'item', 'time_taken', 'distance_travelled', 
         'xPos', 'yPos', 'zPos', 'Pitch', 'Yaw', 'ray.x.', 'ray.y.', 
         'ray.z.', 'ray.distance.', 'distance_goal', 'looking_model', 
         'looking_goal', 'steps', 'correctBlocks', 'incorrectBlocks', 
         'answer')

sc_feat <- sc %>% 
  select(vars) %>% 
  group_by(id) %>%
  summarise(time_taken = max(as.numeric(time_taken), na.rm = T),
            distance_travelled = max(as.numeric(distance_travelled), na.rm = T),
            xPos_range = diff(range(as.numeric(xPos), na.rm = T)), 
            xPos_mean = mean(as.numeric(xPos), na.rm = T), 
            yPos_range = diff(range(as.numeric(yPos), na.rm = T)), 
            yPos_mean = mean(as.numeric(yPos), na.rm = T), 
            zPos_range = diff(range(as.numeric(zPos), na.rm = T)),
            zPos_mean = mean(as.numeric(zPos), na.rm = T), 
            pitch_range = diff(range(as.numeric(Pitch), na.rm = T)),
            pitch_mean = mean(as.numeric(Pitch), na.rm = T), 
            yaw_range = diff(range(as.numeric(Yaw), na.rm = T)),
            yaw_mean = mean(as.numeric(Yaw), na.rm = T),
            ray_x_range = diff(range(as.numeric(ray.x.), na.rm = T)),
            ray_x_mean = mean(as.numeric(ray.x.), na.rm = T),
            ray_y_range = diff(range(as.numeric(ray.y.), na.rm = T)),
            ray_z_mean = mean(as.numeric(ray.z.), na.rm = T),
            ray_dist_range = diff(range(as.numeric(ray.distance.), na.rm = T)),
            ray_dist_mean = mean(as.numeric(ray.distance.), na.rm = T), 
            looking_model = max(as.numeric(looking_model), na.rm = T),
            looking_goal = max(as.numeric(looking_goal), na.rm = T),
            looking_ratio = looking_model/looking_goal)
 
sc_feat %>% head()

sc_feat_item <- sc %>% 
  select(vars) %>% 
  group_by(id, item) %>%
  summarise(time_taken = max(as.numeric(time_taken), na.rm = T),
            distance_travelled = max(as.numeric(distance_travelled), na.rm = T),
            xPos_range = diff(range(as.numeric(xPos), na.rm = T)), 
            xPos_mean = mean(as.numeric(xPos), na.rm = T), 
            yPos_range = diff(range(as.numeric(yPos), na.rm = T)), 
            yPos_mean = mean(as.numeric(yPos), na.rm = T), 
            zPos_range = diff(range(as.numeric(zPos), na.rm = T)),
            zPos_mean = mean(as.numeric(zPos), na.rm = T), 
            pitch_range = diff(range(as.numeric(Pitch), na.rm = T)),
            pitch_mean = mean(as.numeric(Pitch), na.rm = T), 
            yaw_range = diff(range(as.numeric(Yaw), na.rm = T)),
            yaw_mean = mean(as.numeric(Yaw), na.rm = T),
            ray_x_range = diff(range(as.numeric(ray.x.), na.rm = T)),
            ray_x_mean = mean(as.numeric(ray.x.), na.rm = T),
            ray_y_range = diff(range(as.numeric(ray.y.), na.rm = T)),
            ray_z_mean = mean(as.numeric(ray.z.), na.rm = T),
            ray_dist_range = diff(range(as.numeric(ray.distance.), na.rm = T)),
            ray_dist_mean = mean(as.numeric(ray.distance.), na.rm = T), 
            looking_model = max(as.numeric(looking_model), na.rm = T),
            looking_goal = max(as.numeric(looking_goal), na.rm = T),
            looking_ratio = looking_model/looking_goal)


sc_feat_item_wide = sc_feat %>% select(id)
for (i in names(sc_feat_item)[3:length(sc_feat_item)]){
  print(i)
  wide_temp <- sc_feat_item %>% select(id, item, i) %>% spread(key=item, value = i) 
  
  sc_feat_item_wide = sc_feat_item_wide %>% join(wide_temp)
}

sc_feat$id
sc_feat_item_wide


```

```{r}
# correlation table 
sc_feat_all <- sc_feat %>% inner_join(df_all, by="id") 
sc_feat_all %>% select(-id) %>% 
  cor() %>% 
  as.data.frame() %>% 
  rownames_to_column() %>% 
  filter(rowname %in% c("pc", "mr", "sc", "VKMR", "SPM"))

```

# Predict SC Scores From Process Data (Carret )
## Train-test Split
```{r}
library(caret)
library(glmnet)

sc_feat_all_pred_sc <- sc_feat_all %>% select(-c(time_taken, pc, mr, VKMR, SPM, xPos_range, ray_x_range)) %>% column_to_rownames("id")
sc_feat_all_pred_sc %>% head()

in_train <- createDataPartition(y = sc_feat_all_pred_sc$sc,
                                p = 0.7, list = FALSE)

training <- sc_feat_all_pred_sc[ in_train, ]
testing  <- sc_feat_all_pred_sc[-in_train, ]

dim(training)
training

```

## Baseline - Linear Model
```{r}
ctrl <- trainControl(method = "cv", number = 10)

mdl_lm <- caret::train(sc ~ ., data = training, method = "lm", 
             preProcess = c("center", "scale"), 
             trControl = ctrl)

y_pred_lm_cesd <- predict(mdl_lm, newdata = testing)

eval_lm <- defaultSummary(data.frame(obs = testing$sc, pred = y_pred_lm_cesd))
eval_lm
```

## Penalized Linear Model
```{r}
ctrl <- trainControl(method = "cv", number = 10)

Grid <- expand.grid(.lambda = seq(.05, 1, length = 10),
                    .alpha = seq(.05, 1, length = 10))

mdl_enet <- caret::train(sc ~ ., data = training, method = "glmnet",
                  preProcess = c("center", "scale"),
                  trControl = ctrl, tuneGrid = Grid)

y_pred_elnet <- predict(mdl_enet, newdata = testing)

eval_enet <- defaultSummary(data.frame(obs = testing$sc, pred = y_pred_elnet))
eval_enet
```

## Random Forest Regression
```{r}
rf_grid <- data.frame(.mtry = 2:(ncol(training) - 1L))
ctrl <- trainControl(method = "cv", number = 5)

mdl_rf <- caret::train(sc ~ .,
             data = training, method = "rf",
             ntrees = 100, importance = TRUE,
             tuneGrid = rf_grid,
             trControl = ctrl
             )

y_pred_rf <- predict(mdl_rf, newdata = testing)

eval_rf <- defaultSummary(data.frame(obs = testing$sc, pred = y_pred_rf))
eval_rf
```

## Gradient Boostinf
```{r}
bm_grid <- expand.grid(.interaction.depth = seq(1, 7, by = 2),
                        .n.trees = seq(100, 1000, by = 50),
                        .shrinkage = c(0.01, 0.1),
                        .n.minobsinnode = 1:10)

ctrl <- trainControl(method = "cv", number = 5)

mdl_gbm <- caret::train(sc ~ ., data = training, method = "gbm",
             trControl = ctrl, 
             tuneGrid = gbm_grid, 
             verbose = FALSE)

y_pred_gbm <- predict(mdl_gbm, newdata = testing)

eval_gbm <- defaultSummary(data.frame(obs = testing$sc, pred = y_pred_gbm))
eval_gbm
```



# Try MLR
# Predict SC scores
```{r}
library(mlr)

sc_feat_all_pred_sc <- sc_feat_all %>% select(-c(time_taken, pc, mr, VKMR, SPM)) %>% column_to_rownames("id")
sc_feat_all_pred_sc %>% head()

task_sc = TaskRegr$new(id = 'id', backend = sc_feat_all_pred_sc, target = "sc")
task_sc

lrn = lrn("regr.ranger")
rdesc = makeResampleDesc("CV", iters = 5)
res = mlr3::resample(learner = lrn, task = task_sc, resampling = rdesc, measure = list(RMSE))

# train/test split
train_set = sample(task_sc$nrow, 0.8 * task_sc$nrow)
test_set = setdiff(seq_len(task_sc$nrow), train_set)

# train the model
lrn$train(task_sc, row_ids = train_set)

# predict data
prediction = lrn$predict(task_sc, row_ids = test_set)

# calculate performance
prediction$score()

measure = msr("regr.rsq")

# automatic resampling
resampling = rsmp("cv", folds = 10)
rr = resample(task = task_sc, lrn, resampling)
rr$resampling
rr$score(measure)

```

## With nested CV
```{r}
library(paradox)

data("Boston")
Boston

# create task with target sc 
task_sc = TaskRegr$new(id = 'predict_sc', backend = Boston, target = "tax")
task_sc

# get overview of learners
mlr_learners

# create learner with rf algorithm 
lrn_rf = lrn("regr.glmnet")
lrn_rf$param_set

# define inner resampling strategy
rsmp_inner = rsmp("cv", folds=5)

# define evaluation metric
measures = msr("regr.rsq")

# define param grid
param_set = paradox::ParamSet$new(
  params = list(paradox::ParamDbl$new("lambda", lower = 0.001, upper = 0.11)))
param_set

# define stop criterion for parameter search
terminator = term("evals", n_evals = 5)
tuner = tnr("grid_search", resolution = 10)

# define parameter tuning strategy
at_rf = AutoTuner$new(lrn_rf, rsmp_inner, measures = measures, tune_ps = 
  param_set, terminator, tuner = tuner)

# define outer resampling strategy
rsmp_outer = rsmp("cv", folds = 5)
rr_rf = mlr3::resample(task = task_sc, learner = at_rf, resampling = rsmp_outer)

rr_rf$resampling

```


```{r}
library("mlr3tuning")
task = tsk("iris")
task
learner = lrn("classif.rpart")
learner
resampling = rsmp("holdout")
resampling
measures = msr("classif.ce")
measures
param_set = paradox::ParamSet$new(
  params = list(paradox::ParamDbl$new("cp", lower = 0.001, upper = 0.1)))
param_set
terminator = term("evals", n_evals = 5)
terminator
tuner = tnr("grid_search", resolution = 10)
tuner
at = AutoTuner$new(learner, resampling, measures = measures,
  param_set, terminator, tuner = tuner)

resampling_outer = rsmp("cv", folds = 3)
rr = resample(task = task, learner = at, resampling = resampling_outer)
```



## Predict VKMR scores 
```{r}
sc_feat_all_pred_vkmr <- sc_feat_all %>% select(-c(time_taken, pc, mr, sc, SPM)) %>% column_to_rownames("id")
sc_feat_all_pred_vkmr %>% head()

task_vkmr = TaskRegr$new(id = 'id', backend = sc_feat_all_pred_vkmr, target = "VKMR")
task_vkmr

mlr_learners

lrn = lrn("regr.ranger")
rdesc = makeResampleDesc("CV", iters = 5)
res = resample(learner = lrn, task = regr.task, resampling = rdesc, measure = list(RMSE))

# train/test split
train_set = sample(task_vkmr$nrow, 0.8 * task_vkmr$nrow)
test_set = setdiff(seq_len(task_vkmr$nrow), train_set)

# train the model
lrn$train(task_vkmr, row_ids = train_set)

# predict data
prediction = lrn$predict(task_sc, row_ids = test_set)

# calculate performance
prediction$score()

measure = msr("regr.rsq")

# automatic resampling
resampling = rsmp("cv", folds = 10)
rr = resample(task_vkmr, lrn, resampling)
rr$score(measure)
```




## Analyze Narrative Time 
```{r}

pc_narrative <- pc %>% 
  select(id, narrative_time) %>% 
  rename(pc_narrative_time = narrative_time) %>% 
  group_by(id) %>% 
  summarise(pc_narrative_time = max(pc_narrative_time))

mr_narrative <- mr %>% 
  select(id, narrative_time) %>% 
  rename(mr_narrative_time = narrative_time) %>% 
  group_by(id) %>% 
  summarise(mr_narrative_time = max(mr_narrative_time))

sc_narrative <- sc %>% 
  select(id, narrative_time) %>% 
  rename(sc_narrative_time = narrative_time) %>% 
  group_by(id) %>% 
  summarise(sc_narrative_time = max(sc_narrative_time))

df_narrative = join_all(list(df_all, pc_narrative, mr_narrative, sc_narrative), by = 'id')

df_narrative <- df_narrative %>% mutate(sum_narrative_time  = pc_narrative_time+mr_narrative_time+sc_narrative_time)

cor(df_narrative)

```