# Set working Directory:
setwd("C:/perbb/Chicago2014")

# Import data:
sensintro <- read.table("introexample.csv",header=TRUE,sep=",",dec=".")

# Make relevant variables into factors:
sensintro$assessor=factor(sensintro$assessor)
## # Check your data:
 head(sensintro)
 summary(sensintro)
# Check your data:
head(sensintro)
summary(sensintro)
# Make a data file with only the first 10 assessors
sensintro10=subset(sensintro,as.integer(assessor)<11)

# Fixed effects analysis:
anova(lm(sweet~assessor+product,data=sensintro))

# Same as paired t-test on the 10 assessors with no missings:
with(sensintro10, t.test(sweet~product,paired=TRUE))
# The mixed model version of it:
# We first just load the lme4 package 
# If first time: first install the package
library(lme4)

# Analysing the 10 first assessors:
lmer0 <- lmer(sweet~product+(1|assessor), data=sensintro10)
summary(lmer0)
anova(lmer0)


# Analysing all 20 assessors:
lmer1 <- lmer(sweet~product+(1|assessor), data=sensintro)
summary(lmer1)
anova(lmer1)


# We next load the lmerTest package, 
# which actually also loads the lme4-package 
# So in the future you don't need the step above
# If first time: first install the package
library(lmerTest)

# Now check the anova-output:
lmer1 <- lmer(sweet~product+(1|assessor), data=sensintro)
anova(lmer1)
# Or to get (almost) everything - 
# use the step-function of lmerTest:
# We use it now WITHOUT the automated model reduction feature
s <- step(lmer1,reduce.fixed = FALSE, reduce.random = FALSE)
s
# To check what the R-names of the various output is:
# str(s)
s$diffs.lsmeans.table
plot(s)
# From Kuznetsova et al (2014): 7.1.2. Automated analysis of TVbo data set in R
# Running the largest possible model:
modelTVbo <- lmer(Colourbalance~TVset*Picture+(1|Assessor) + (1 | Assessor:TVset)
                  + (1 | Assessor:Picture) + (1 | Assessor:TVset:Picture)+ (1 | Repeat) 
                  + (1 | Repeat:TVset) + (1 | Repeat:Picture)+ (1 | Repeat:TVset:Picture)
                  ,data=TVbo)

# Running the step-function of lmerTest with default options:
stepTVbo <- step(modelTVbo)

# Look at (all) the results
stepTVbo
# Look at the (default) plots
plot(stepTVbo)
## # Showing ALL the (default)options for the step function:
## stepTVbo=step(modelTVbo, ddf="Satterthwaite", type=3, alpha.random = 0.1,
##               alpha.fixed = 0.05, reduce.fixed = TRUE, reduce.random = TRUE,
##               lsmeans.calc=TRUE, difflsmeans.calc=TRUE, test.effs=NULL,
##               method.grad="simple")
## # Check the details: (choose the help info regarding the lmerTest version of step )
## ?step
# Now the ham data:
data(ham)
modelHam<-lmer(Informed.liking ~ Product*Information*Gender*Age+
(1|Consumer) +(1|Product:Consumer) + (1|Information:Consumer), data=ham)

stepHam<-step(modelHam)

# Final model
stepHam$model

# Overview of random effects
stepHam$rand.table

# Overview of fixed effects
stepHam$anova.table

# Differences of LSMEANS for final model 
stepHam$diffs.lsmeans.table
# Look at the (default) plots
plot(stepHam)
# And finally the carrots data:
modelCarrots <- lmer(Preference~sens2*sens1*Homesize*Age+(1 | product) + (1 + sens1 + sens2 | Consumer), data=carrots)

stepCarrots=step(modelCarrots)

stepCarrots
# Look at the (default) plots
plot(stepCarrots)