### R code from vignette source 'vignettes/PhenStat/inst/doc/PhenStat.Rnw'

###################################################
### code chunk number 1: R.hide
###################################################
library(PhenStat)

dataset1 <- system.file("extdata", "test1.csv", package="PhenStat")

dataset2 <- system.file("extdata", "test1.txt", package="PhenStat")


###################################################
### code chunk number 2: R.hide
###################################################
# Default behaviour with messages
library(PhenStat)
dataset1 <- system.file("extdata", "test1.csv", package="PhenStat")
test <- PhenList(dataset=read.csv(dataset1),
        testGenotype="Sparc/Sparc")

# Out-messages are switched off 
test <- PhenList(dataset=read.csv(dataset1),
        testGenotype="Sparc/Sparc",
        outputMessages=FALSE)


###################################################
### code chunk number 3: R.hide
###################################################
library(PhenStat)
dataset1 <- system.file("extdata", "test3.csv", package="PhenStat")

test <- PhenList(dataset=read.csv(dataset1), 
        dataset.clean=TRUE, 
        dataset.values.female=1, 
        dataset.values.male=2, 
        testGenotype="Mysm1/+")


###################################################
### code chunk number 4: R.hide
###################################################
library(PhenStat)
dataset1 <- system.file("extdata", "test3.csv", package="PhenStat")

test <- PhenList(dataset=read.csv(dataset1), 
        dataset.clean=TRUE, 
        dataset.values.female=1, 
        dataset.values.male=2, 
        testGenotype="Mysm1/+")

test$dataset

test$dataset.stat


###################################################
### code chunk number 5: R.hide
###################################################
library(PhenStat)
dataset2 <- system.file("extdata", "test2.csv", package="PhenStat")

test2 <- PhenList(dataset=read.csv(dataset2),
        testGenotype="Arid4a/Arid4a",
        dataset.colname.weight="Weight.Value")

test2$testGenotype

test2$refGenotype

test2$dataset.colname.weight


###################################################
### code chunk number 6: R.hide
###################################################
library(PhenStat)
dataset1 <- system.file("extdata", "test1.csv", package="PhenStat")

test <- PhenList(dataset=read.csv(dataset1),
        testGenotype="Sparc/Sparc", 
        outputMessages=FALSE)

# Default behaviour
result <- testDataset(test,
        depVariable="Bone.Area", 
        equation="withoutWeight")

# Perform each step of the MM framework separatly
result <- testDataset(test,
        depVariable="Bone.Area", 
        equation="withoutWeight",callAll=FALSE)

# Estimated model effects
result$model.effect.batch

result$model.effect.variance

result$model.effect.weight

result$model.effect.sex

result$model.effect.interaction

result$numberSexes

# Change the effect values: interaction effect will stay in the model
result <- testDataset(test,
        depVariable="Bone.Area", 
        equation="withoutWeight",
        keepList=c(TRUE,TRUE,FALSE,TRUE,TRUE),
        callAll=FALSE)

result <- finalModel(result)

summaryOutput(result)


###################################################
### code chunk number 7: R.hide
###################################################
testFinalModel(result)

classificationTag(result)


###################################################
### code chunk number 8: R.hide
###################################################
file <- system.file("extdata", "test7_TFE.csv", package="PhenStat")
test <- PhenList(dataset=read.csv(file),
        testGenotype="het",
        refGenotype = "WT",
        dataset.colname.sex="sex",
        dataset.colname.genotype="Genotype",
        dataset.values.female="f",
        dataset.values.male= "m",
        dataset.colname.weight="body.weight",
        dataset.colname.batch="Date_of_procedure_start")

# TFDataset function creates cleaned dataset - concurrent controls dataset
test_TF <- TFDataset(test,depVariable="Cholesterol")

# TF method is called
result  <- testDataset(test_TF,
        depVariable="Cholesterol",
        method="TF")
summaryOutput(result)


###################################################
### code chunk number 9: R.hide
###################################################
library(PhenStat)
file <- system.file("extdata", "test1.csv", package="PhenStat")
test <- PhenList(dataset=read.csv(file),
        testGenotype="Sparc/Sparc")

# RR method is called
result <- testDataset(test,
        depVariable="Lean.Mass",
        method="RR")
summaryOutput(result)


###################################################
### code chunk number 10: R.hide
###################################################
library(PhenStat)
dataset_cat <- system.file("extdata", "test_categorical.csv", 
                           package="PhenStat")

test_cat <- PhenList(read.csv(dataset_cat),testGenotype="Aff3/Aff3")

result_cat <- testDataset(test_cat,
        depVariable="Thoracic.Processes",
        method="FE")

result_cat$depVariable

result_cat$method

result_cat$numberSexes

# Chi squared table for all data
result_cat$model.output$count_matrix_all

# Chi squared table for males only records
result_cat$model.output$count_matrix_male

# Percentage matrix for all data
result_cat$model.output$percentage_matrix_all

# Percentage matrix for females only records
result_cat$model.output$percentage_matrix_female

# Matrix statistics for all data
result_cat$model.output$stat_all

# Matrix statistics for males only records
result_cat$model.output$stat_male

# Effect size for all data
result_cat$model.output$ES

# Effect size for females only records
result_cat$model.output$ES_female

# Fisher Exact Test results for all data
result_cat$model.output$all

# p-value for all data
result_cat$model.output$all$p.value


###################################################
### code chunk number 11: R.hide
###################################################
library(PhenStat)
dataset1 <- system.file("extdata", "test1.csv", package="PhenStat")

# MM framework
test <- PhenList(dataset=read.csv(dataset1),
        testGenotype="Sparc/Sparc",outputMessages=FALSE)

result <- testDataset(test,
        depVariable="Lean.Mass",
        outputMessages=FALSE)

summaryOutput(result)


###################################################
### code chunk number 12: R.hide
###################################################
library(PhenStat)
dataset_cat <- system.file("extdata", "test_categorical.csv", 
                           package="PhenStat")

test2 <- PhenList(dataset=read.csv(dataset_cat),
        testGenotype="Aff3/Aff3",outputMessages=FALSE)

result2 <- testDataset(test2,
        depVariable="Thoracic.Processes",
        method="FE",outputMessages=FALSE)  

summaryOutput(result2)


###################################################
### code chunk number 13: R.hide
###################################################
library(PhenStat)
dataset_cat <- system.file("extdata", "test_categorical.csv", 
                           package="PhenStat")

test_cat <- PhenList(dataset=read.csv(dataset_cat),
        testGenotype="Aff3/Aff3",outputMessages=FALSE)

result_cat <- testDataset(test_cat,
        depVariable="Thoracic.Processes",
        method="FE",outputMessages=FALSE)  

vectorOutput(result_cat)


###################################################
### code chunk number 14: R.hide
###################################################
library(PhenStat)
dataset_cat <- system.file("extdata", "test_categorical.csv", 
                           package="PhenStat")

test_cat <- PhenList(dataset=read.csv(dataset_cat),
        testGenotype="Aff3/Aff3",outputMessages=FALSE)

result_cat <- testDataset(test_cat,
        depVariable="Thoracic.Processes",
        method="FE",outputMessages=FALSE) 

vectorOutputMatrices(result_cat)


###################################################
### code chunk number 15: R.hide
###################################################
library(PhenStat)
dataset_cat <- system.file("extdata", "test_categorical.csv", 
                           package="PhenStat")

test_cat <- PhenList(dataset=read.csv(dataset_cat),
        testGenotype="Aff3/Aff3",outputMessages=FALSE)

result_cat <- testDataset(test_cat,
        depVariable="Thoracic.Processes",
        method="FE",outputMessages=FALSE) 

categoricalBarplot(result_cat)


###################################################
### code chunk number 16: R.hide
###################################################
library(PhenStat)
dataset1 <- system.file("extdata", "test1.csv", package="PhenStat")

# MM framework
test <- PhenList(dataset=read.csv(dataset1),
        testGenotype="Sparc/Sparc",outputMessages=FALSE)

result <- testDataset(test,
        depVariable="Lean.Mass",
        outputMessages=FALSE)

boxplotSexGenotype(test,
        depVariable="Lean.Mass",
        graphingName="Lean Mass")

scatterplotSexGenotypeBatch(test,
        depVariable="Lean.Mass",
        graphingName="Lean Mass")

scatterplotGenotypeWeight(test,
        depVariable="Bone.Mineral.Content",
        graphingName="BMC")


###################################################
### code chunk number 17: R.hide
###################################################
library(PhenStat)
dataset1 <- system.file("extdata", "test1.csv", package="PhenStat")

# MM framework
test <- PhenList(dataset=read.csv(dataset1),
        testGenotype="Sparc/Sparc",outputMessages=FALSE)

result <- testDataset(test,
        depVariable="Lean.Mass",
        outputMessages=FALSE)

qqplotGenotype(result)

qqplotRandomEffects(result)

qqplotRotatedResiduals(result)

plotResidualPredicted(result)

boxplotResidualBatch(result)