## ----echo=FALSE---------------------------------------------------------------
library(DACF)
library(qgraph)
r.1=rnorm(1000,20,5)
r.2=rnorm(1000,30,5)
r.3=rnorm(1000,10,10)
ca.true=matrix(c(r.1,r.2,r.3,rep(c(2,1,3),each=1000)),ncol=2)
r.2.cf=induce.cfe(0,.3,r.2)
r.3.cf=induce.cfe(.3,0,r.3)
ca.cf=matrix(c(r.1,r.2.cf,r.3.cf,rep(c(2,1,3),each=1000)),ncol=2)
colnames(ca.cf)=c('score','group')

## -----------------------------------------------------------------------------
# group sample mean
aggregate(ca.true[,1],mean,by=list(ca.true[,2]))
# group sample variance
aggregate(ca.true[,1],var,by=list(ca.true[,2]))

## -----------------------------------------------------------------------------
# group sample mean
aggregate(ca.cf[,1],mean,by=list(ca.cf[,2]))
# group sample variance
aggregate(ca.cf[,1],var,by=list(ca.cf[,2]))

## -----------------------------------------------------------------------------
# younger-aged group
young=ca.cf[ca.cf[,2]==1,1]
rec.mean.var(young, flr=min(young), ceil = max(young)) # true mean and variance are 30 and 25

# the estimated floor and ceiling percentages and the recovered mean and variance estimates are displayed above

# older-aged group
old=ca.cf[ca.cf[,2]==3,1]
rec.mean.var(old, flr=min(old), ceil = max(old)) # true mean and variance are 10 and 100

# the estimated floor and ceiling percentages and the recovered mean and variance estimates are displayed above

## -----------------------------------------------------------------------------
ca.cf.df = data.frame(ca.cf)

# ANOVA
lw.f.star(ca.cf.df, score ~ group, flr = min(ca.cf.df$score), ceil = max(ca.cf.df$score), "a")
lw.f.star(ca.cf.df, score ~ group, flr = min(ca.cf.df$score), ceil = max(ca.cf.df$score), "b")

# multiple comparisons after ANOVA
aov.cf=aov(score ~ group, data=ca.cf.df)
cf_correction(aov.cf, tests="all", df.adjustment="trunc",
              gh.correction="no_gh", alpha=.05,
              flr=min(ca.cf.df$score), ceil=max(ca.cf.df$score))

# t-test
lw.t.test(young, old, flr1 = min(young), ceil1 = max(young), flr2 = min(old), ceil2 = max(old), "a")
lw.t.test(young, old, flr1 = min(young), ceil1 = max(young), flr2 = min(old), ceil2 = max(old), "b")

## -----------------------------------------------------------------------------
# Simulate healthy data for two groups
x.1=rnorm(300,2,4)
x.2=rnorm(300,3,5)
# check mean and variance for simulated healthy data
mean(x.1);var(x.1)
mean(x.2);var(x.2)
# induce ceiling effects of 20% in group 1
x.1.cf=induce.cfe(.2,0,x.1)
# induce floor effects of 10% in group 2
x.2.cf=induce.cfe(0,.1,x.2)
# recover the mean and variance for ceiling/floor data
rec.mean.var(x.1.cf, flr = min(x.1.cf), ceil = max(x.1.cf))
rec.mean.var(x.2.cf, flr = min(x.2.cf), ceil = max(x.2.cf))
# conduct a t test on healthy data
t.test(x.1,x.2)
t.test(x.1.cf,x.2.cf)
# conduct an adjusted t test on ceiling/floor data
lw.t.test(x.1.cf, x.2.cf, flr1 = min(x.1.cf), ceil1 = max(x.1.cf), flr2 = min(x.2.cf), ceil2 = max(x.2.cf), "a")
lw.t.test(x.1.cf, x.2.cf, flr1 = min(x.1.cf), ceil1 = max(x.1.cf), flr2 = min(x.2.cf), ceil2 = max(x.2.cf), "b")
# generate a dataframe for ANOVA demo
testdat=threeganova.sim(10000,.0625,1)
# induce ceiling/floor effects in the data
testdat.cf=testdat
testdat.cf[testdat.cf$group==2,]$y=induce.cfe(.2,0,testdat.cf[testdat.cf$group==2,]$y)
# conduct an adjusted F star test on ceiling/floor data
lw.f.star(testdat.cf, y ~ group, flr = min(testdat.cf$y), ceil = max(testdat.cf$y), "a")
lw.f.star(testdat.cf, y ~ group, flr = min(testdat.cf$y), ceil = max(testdat.cf$y), "b")
# conduct post-hoc comparisons with Games-Howell/Welch correction
aov.testdat.cf=aov(y ~ group, data=testdat.cf)
cf_correction(aov.testdat.cf, tests="all", df.adjustment="trunc",
              gh.correction="yes_gh", alpha=.05,
              flr=min(testdat.cf$y), ceil=max(testdat.cf$y))

## ----eval=T-------------------------------------------------------------------

set.seed(123)
Sigma=matrix(.30,6,6)
diag(Sigma)=1
y.true=MASS::mvrnorm(400,rep(0,6),Sigma)
colnames(y.true)=paste0("V",1:6)

y.cf=pmin(pmax(y.true,-1.5),1.5)

fit=ggm_cfe(y.cf, floor=-1.5, ceiling=1.5, method="EBICglasso")
print(fit)
summary(fit)
plot(fit)

## ----eval=T-------------------------------------------------------------------
Sigma.corrected=fit$Sigma_corrected

qgraph(Sigma.corrected,
       graph="glasso",
       sampleSize=fit$n,
       gamma=.5,
       layout="spring")

qgraph(Sigma.corrected,
       graph="pcor",
       threshold="fdr",
       sampleSize=fit$n,
       layout="spring")

## ----eval=T-------------------------------------------------------------------
fit.all=ggm_cfe(y.cf, floor=-1.5, ceiling=1.5, method="all")

sapply(fit.all$networks_corrected, function(x) sum(x[upper.tri(x)]!=0))
sapply(fit.all$networks_naive, function(x) sum(x[upper.tri(x)]!=0))