---
title: "Data Analyses with Ceiling/Floor data"
author: "Qimin Liu"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Data Analyses with Ceiling/Floor data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Summary 

Ceiling and floor effects are common in data. Ceiling or floor effects occur when the tests or scales are relatively easy or difficult such that substantial proportions of individuals obtain either maximum or minimum scores and that the true extent of their abilities cannot be determined. 

Ceiling and floor effects, subsequently, causes problems in data analysis. For example, ceiling or floor effects alone would induce, respectively, attenuation or inflation in mean estimates. And both ceiling and floor effects would result in attenuation in variance estimates. This imposes challenges in mean and variance based data analytic methods.

This package implements methods to deal with challenges associated with ceiling/floor effects in the data using parametric methods that assume normality for the true scores. The current version is capable of mean and variance recovery given data with ceiling/floor effects and of mean comparison tests such as t-test and ANOVA for data with ceiling/floor effects.

# Helper functions

The package contains a helper function `threeganova.sim` that would generate a three-group anova data with a standard normal control group and positive/negative treatment groups of effect with same magnitudes. In addition, one can specify the standard deviation in positive treatment group. To see the specifics of the function, user can enter `?threeganova.sim` in the R console.

Another helper function included in the package is `induce.cfe` where the user can manually induce ceiling and floor effects to healthy data. To see the specifics of the function, user can enter `?induce.cfe` in the R console.

Moreover, the function `F.star.test` allows user to conduct a Brown-Forsythe F star test. This is a variant of the commonly used F test. F star test is robust against violations of homogeneity of variance (HOV) assumption for the F test.

# Functions for data analyses

The package includes several functions that can facilitate the user to conduct data analyses for data with ceiling/floor effects. `rec.mean.var` estimates the true mean and variance of the data with ceiling/floor effects. That is, as mentioned in the summary, the observed mean and variance of data with ceiling/floor effects are often biased. Thus, `rec.mean.var` aims to help the user to recover the mean and variance of the data were ceiling/floor effects absent. `lw.t.test` conducts a t test that adjusts for ceiling/floor effects in the data. As `lw.t.test` also uses Welch's t test, the adjusted t test is robust against HOV violation. `lw.f.star` conducts a F star test for one-way ANOVA that adjusts for ceiling/floor effects in the data. `lw.f.star` is also robust against HOV violation. For both `lw.f.star` and `lw.t.test`: method `a` is a liberal approach that yields accurate effect size estimates but has mildly inflated type I error rates, `b` is a conservative approach with well-controlled type I error rates that have good, but less accurate than `a`, effect estimates.

As a quick guide:

- Use `rec.mean.var` to recover a mean and variance for one variable.
- Use `lw.t.test` to compare two groups.
- Use `lw.f.star` to conduct a one-way ANOVA.
- Use `cf_correction` for post-hoc pairwise comparisons after ANOVA.
- Use `ggm_cfe` when the variables are analyzed as a psychological network or Gaussian graphical model.

When the omnibus ANOVA is followed by pairwise group comparisons, `cf_correction` extends Tukey's HSD and Bonferroni procedures to account for ceiling/floor effects. The function can use either a dataframe with group IDs and scores or an `aov` object. For dataframe input, the first column should contain group IDs and the second column should contain the outcome values. For clarity, the examples below use `aov` objects. The function returns adjusted group means, pairwise mean differences, confidence intervals, p values, and effect sizes.

For multivariate network data, `ggm_cfe` estimates a Gaussian graphical model after correcting the covariance matrix for ceiling/floor effects. Rows should be participants or observations and columns should be variables, items, or symptoms. The function first recovers marginal means and variances, then estimates pairwise latent correlations by method of moments, and finally supplies the corrected covariance matrix to `qgraph` for EBICglasso, FDR, or Bonferroni edge selection.

The `flr` and `ceil` arguments in `rec.mean.var`, `lw.t.test`, `lw.f.star`, and `cf_correction` should be set to the known minimum and maximum of the measurement scale (e.g., 0 and 100 for a percentage scale). The same idea applies to the `floor` and `ceiling` arguments in `ggm_cfe`. When the scale bounds are unknown, `min()` and `max()` of the observed data can be used as reasonable approximations.

# Example 1: an Aging Example on t-test, ANOVA, and Tukey's HSD

Imagine a scenario where we wish to test the difference in cognitive ability for people of different age groups. In this toy example, we have 1000 participants for three age groups, the younger-aged group has true mean and variance of respectively 30 and 25, the middle-aged group 20 and 25 and the older-aged group 10 and 100. The higher the score, the higher the cognitive ability. We can check the mean and variance of the true mean and variance on the data composed of true scores, `ca.true`. 

```{r 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')
```

```{r}
# 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]))
```

Now consider the fact that a substantial proportion of the younger-aged group may score maximum at the cognitive ability test and a substantial proportion of the older-aged group may score minimum. Let both the ceiling and the floor proportions be 15%, we have the dataset `ca.cf`.

```{r}
# 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]))
```

We can see that both the mean and the variance estimates from the younger-aged and the older-aged groups are biased. The function `rec.mean.var` can help recover the mean and variance. In the example of the younger-aged group, we first select all the scores of the younger-aged group and name it as a new variable `young` and then use our function `rec.mean.var` to recover the mean and variance. We can do the same for the older-aged group.

```{r}
# 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
```

Now we wish to conduct an ANOVA in the data with floor and ceiling effects. We can use the function `lw.f.star`. We can also conduct a t-test between the older-aged and the younger-aged group by using the function `lw.t.test`. Both methods `a` and `b` are used for the illustration purposes.

```{r}
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")
```
The ANOVA output tests whether there are any group differences. The `cf_correction` output then shows the adjusted pairwise comparisons. The columns `Comparison_i` and `Comparison_j` identify the two groups being compared, `diff_in_means` gives the adjusted mean difference, and the Tukey and Bonferroni columns give the corresponding confidence intervals and p values. Both the ANOVA and the t-tests returned significant results. 

# Example 2: A Simulated Data Example for t-test, ANOVA, and Games-Howell Pairwise Comparison

The following example provides an overview of the helper functions in the package that can aid in simulations and further demonstrates data analytic functions in the package.

```{r}
# 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))
```

The final call requests both Tukey and Bonferroni pairwise comparisons and includes the Games-Howell/Welch correction for unequal variances. When `gh.correction="yes_gh"`, the output includes additional Games-Howell and Welch columns, which are useful when the group variances are not assumed to be equal.

# Example 3: Network Analysis with Ceiling/Floor Data

The function `ggm_cfe` can be used when the variables are items or symptoms that will be analyzed as a Gaussian graphical model. The input should be a numeric matrix or dataframe with one row per person and one column per variable. The `floor` and `ceiling` arguments should be set to the lower and upper limits of the scale. In the following example, we simulate six correlated variables, impose floor and ceiling effects, and estimate the corrected network.

```{r 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)
```

The output first reports the proportion of observations at the floor and ceiling for each variable and the size of the nearest positive definite correction. The object `fit$network_corrected` contains the corrected partial correlation network, and `fit$network_naive` contains the corresponding network from the uncorrected observed covariance matrix. The network plot shows the corrected network alongside the naive network that ignores ceiling/floor effects.

Users can also extract the corrected covariance matrix and pass it directly to `qgraph`. This is useful when the user wishes to use familiar `qgraph` syntax for network estimation, plotting, and centrality. The corrected covariance matrix is stored in `fit$Sigma_corrected`.

```{r 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")
```

The first `qgraph` call estimates an EBICglasso network from the corrected covariance matrix. The second call estimates a partial correlation network and retains edges using FDR correction. Other `qgraph` options, such as node labels, layouts, and centrality functions, can then be used in the usual way.

The same corrected covariance matrix can also be used with significance-based edge selection inside `ggm_cfe`. Setting `method="all"` returns EBICglasso, FDR, and Bonferroni networks in the same object.

```{r 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))
```

The first line gives the number of retained edges after the ceiling/floor correction. The second line gives the corresponding edge counts from the naive analysis.
