#####################
# This R script can be used to fit linear mixed models for the different
# dependent measures. You only need to change the working directory to the 
# path you extracted the data to.
# Further, two additional packages are required:
# lme4 (performs the actual fitting)
# lmerTest (allows convenient statistic comparisons)
# You can install these packages via the shell (internet connection required):
# install.packages("lme4", dependencies = TRUE)
# install.packages("lmerTest", dependencies = TRUE)
#
# The data contains the three dependent measures: palm drift, angular disparity
# and positional discrepancy. The analysis is separated for the three measures.
# The procedure is always the same:
# - find the null model with the maximal, converging random effect structure 
# - compare increasingly complex models to find the relevant fixed effects
# - a confirmatory ANOVA performed with the lmerTest package on the full model
#####################

library(lme4)
library(lmerTest)

# enter the working directory here
setwd("")
# read in the data
data <- read.table("LiS_LME_Data.dat", header = TRUE)
# just to check if everything is there
summary(data)

##################### Palm Drift

########### Null model selection

pd.model.nullfull = lmer(PalmDrift ~ 1 + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE) # no convergence

# single effect null models
pd.model.nullvis = lmer(PalmDrift ~ 1 + (1+HandVisibility|ParticipantID), data = data, REML = FALSE)
pd.model.nullref = lmer(PalmDrift ~ 1 + (1+Reference|ParticipantID), data = data, REML = FALSE)
pd.model.nulloff = lmer(PalmDrift ~ 1 + (1+OffsetType|ParticipantID), data = data, REML = FALSE) # no convergence

# two effect null models
pd.model.nullvisref = lmer(PalmDrift ~ 1 + (1+HandVisibility|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE) # no convergence
pd.model.nullrefoff = lmer(PalmDrift ~ 1 + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
pd.model.nullvisoff = lmer(PalmDrift ~ 1 + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE) # best likelihood

# likelihoodratio tests
anova(pd.model.nullvis, pd.model.nullrefoff)
anova(pd.model.nullref, pd.model.nullrefoff)

anova(pd.model.nullvis, pd.model.nullvisoff)
anova(pd.model.nullref, pd.model.nullvisoff)

anova(pd.model.nullrefoff, pd.model.nullvisoff)

# model applying random slopes / random intercepts for hand visibility and offset type provides best fit and has the highest loglikelihood of the considered two factor models
plot(fitted(pd.model.nullvisoff), residuals(pd.model.nullvisoff))
coef(pd.model.nullvisoff)
summary(pd.model.nullvisoff)

########### Fixed effect selection
# single effect models
pd.model.vis = lmer(PalmDrift ~ HandVisibility + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE) # best model
pd.model.off = lmer(PalmDrift ~ OffsetType + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
pd.model.ref = lmer(PalmDrift ~ Reference + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)

anova(pd.model.nullvisoff, pd.model.vis) # significant improvement
anova(pd.model.nullvisoff, pd.model.ref)
anova(pd.model.nullvisoff, pd.model.off)
# complex models
pd.model.visoff = lmer(PalmDrift ~ HandVisibility + OffsetType + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
pd.model.visref = lmer(PalmDrift ~ HandVisibility + Reference + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
pd.model.full   = lmer(PalmDrift ~ HandVisibility * Reference * OffsetType + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)

anova(pd.model.vis, pd.model.visoff)
anova(pd.model.vis, pd.model.visref) # nearly significant improvement
anova(pd.model.vis, pd.model.full)

# single effect model using hand visibility as fixed effect provides best and most parsimonious fit

summary(pd.model.vis)
anova(pd.model.vis)
rand(pd.model.vis)

########### Full model ANOVA

summary(pd.model.full)
anova(pd.model.full)
rand(pd.model.full)

##################### Angular Disparity

########### Null model selection

ad.model.nullfull = lmer(AngularDisparity ~ 1 + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE) # converges

# the full null model applying random slopes and intercept for all factors converges and can be used as baseline

########### Fixed effect selection
# single effect models
ad.model.vis = lmer(AngularDisparity ~ HandVisibility + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)
ad.model.off = lmer(AngularDisparity ~ OffsetType + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)
ad.model.ref = lmer(AngularDisparity ~ Reference + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)

anova(ad.model.nullfull, ad.model.vis) # significant improvement
anova(ad.model.nullfull, ad.model.off)
anova(ad.model.nullfull, ad.model.ref) # significant improvement

# best models are those assuming a fixed effect for either reference or hand visibility

# more complex models
ad.model.visref    = lmer(AngularDisparity ~ HandVisibility + Reference + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)
ad.model.visrefia  = lmer(AngularDisparity ~ HandVisibility * Reference + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)
ad.model.visrefoff = lmer(AngularDisparity ~ HandVisibility + Reference + OffsetType + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE) # no convergence
ad.model.full      = lmer(AngularDisparity ~ HandVisibility * Reference * OffsetType + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)


data.model.visref.m = lmer(AngularDisparity ~ HandVisibility + Reference + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE)
summary(data.model.visref.m)
anova(data.model.visref.m)
rand(data.model.visref.m)

anova(ad.model.vis, ad.model.visref) # significant improvement
anova(ad.model.ref, ad.model.visref) # significant improvement
anova(ad.model.visref, ad.model.visrefia) # nearly significant improvement
anova(ad.model.visref, ad.model.full)

# two effect model using hand visibility and reference as fixed effect provides best and most parsimonious fit

summary(ad.model.visref)
anova(ad.model.visref)
rand(ad.model.visref)

########### Full model ANOVA

summary(ad.model.full)
anova(ad.model.full)
rand(ad.model.full)

##################### Positional Discrepancy

########### Null model selection

dis.model.nullfull = lmer(PositionDiscrepancy ~ 1 + (1+HandVisibility|ParticipantID) + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE) # no convergence

# single effect null models
dis.model.nullvis = lmer(PositionDiscrepancy ~ 1 + (1+HandVisibility|ParticipantID), data = data, REML = FALSE)
dis.model.nullref = lmer(PositionDiscrepancy ~ 1 + (1+Reference|ParticipantID), data = data, REML = FALSE)
dis.model.nulloff = lmer(PositionDiscrepancy ~ 1 + (1+OffsetType|ParticipantID), data = data, REML = FALSE) # no convergence

# two effect null models
dis.model.nullvisref = lmer(PositionDiscrepancy ~ 1 + (1+HandVisibility|ParticipantID) + (1+Reference|ParticipantID), data = data, REML = FALSE) # no convergence
dis.model.nullrefoff = lmer(PositionDiscrepancy ~ 1 + (1+Reference|ParticipantID)      + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
dis.model.nullvisoff = lmer(PositionDiscrepancy ~ 1 + (1+HandVisibility|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE) # no convergence

# likelihoodratio tests

anova(dis.model.nullvis, dis.model.nullrefoff)
anova(dis.model.nullref, dis.model.nullrefoff)

# model applying random slopes / random intercepts for hand reference and offset type provides best fit of the considered null models
plot(fitted(dis.model.nullrefoff), residuals(dis.model.nullrefoff))
coef(dis.model.nullrefoff)
summary(dis.model.nullrefoff)


########### Fixed effect selection
# single effect models
dis.model.vis = lmer(PositionDiscrepancy ~ HandVisibility + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
dis.model.off = lmer(PositionDiscrepancy ~ OffsetType     + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
dis.model.ref = lmer(PositionDiscrepancy ~ Reference      + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)

anova(dis.model.nullrefoff, dis.model.vis) # significant improvement
anova(dis.model.nullrefoff, dis.model.off)
anova(dis.model.nullrefoff, dis.model.ref) # significant improvement

# models applying visibility or reference as fixed factors improve the fit significantly
# two effect model
dis.model.visref = lmer(PositionDiscrepancy ~ HandVisibility + Reference  + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)

anova(dis.model.visref, dis.model.vis) # significant improvement
anova(dis.model.visref, dis.model.ref) # significant improvement

# two effect model improves fit significantly compared to both of the single effect models
# complex models

dis.model.visoffiaref  = lmer(PositionDiscrepancy ~ HandVisibility * OffsetType + Reference + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
dis.model.visrefoffia  = lmer(PositionDiscrepancy ~ HandVisibility + Reference * OffsetType + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
dis.model.full         = lmer(PositionDiscrepancy ~ HandVisibility * OffsetType * Reference + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)

data.model.visoffiaref.m = lmer(PositionDiscrepancy ~ HandVisibility * OffsetType + Reference + (1+Reference|ParticipantID) + (1+OffsetType|ParticipantID), data = data, REML = FALSE)
summary(data.model.visoffiaref.m)
anova(data.model.visoffiaref.m)
rand(data.model.visoffiaref.m)


anova(dis.model.visref, dis.model.visoffiaref) # significant improvement
anova(dis.model.visref, dis.model.visrefoffia)
anova(dis.model.visoffiaref, dis.model.full)

# model using hand visibility and reference as fixed effect, as well as the interaction between visibility and offset provides best and most parsimonious fit; the main effect of offset type is not significant however
summary(dis.model.visoffiaref)
anova(dis.model.visoffiaref)
rand(dis.model.visoffiaref)

########### Full model ANOVA

summary(dis.model.full)
anova(dis.model.full)
rand(dis.model.full)