Data Analyses with Ceiling/Floor data

Qimin Liu

2026-07-16

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:

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.

# group sample mean
aggregate(ca.true[,1],mean,by=list(ca.true[,2]))
##   Group.1         x
## 1       1 30.087152
## 2       2 20.026629
## 3       3  9.254644
# group sample variance
aggregate(ca.true[,1],var,by=list(ca.true[,2]))
##   Group.1         x
## 1       1  25.16958
## 2       2  25.73606
## 3       3 106.66159

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.

# group sample mean
aggregate(ca.cf[,1],mean,by=list(ca.cf[,2]))
##   Group.1        x
## 1       1 29.12765
## 2       2 20.02663
## 3       3 11.21187
# group sample variance
aggregate(ca.cf[,1],var,by=list(ca.cf[,2]))
##   Group.1        x
## 1       1 14.18011
## 2       2 25.73606
## 3       3 59.81278

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.

# 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
## $ceiling.percentage
## [1] 0.298
## 
## $floor.percentage
## [1] 0
## 
## $est.mean
## [1] 30.06951
## 
## $est.var
## [1] 24.52914
# 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
## $ceiling.percentage
## [1] 0
## 
## $floor.percentage
## [1] 0.307
## 
## $est.mean
## [1] 11.30189
## 
## $est.var
## [1] 38.30317
# 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.

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")
## $statistic
## [1] 2450.989
## 
## $p.value
## [1] 0
## 
## $est.f.squared
## [1] 1.633993
lw.f.star(ca.cf.df, score ~ group, flr = min(ca.cf.df$score), ceil = max(ca.cf.df$score), "b")
## $statistic
## [1] 2450.989
## 
## $p.value
## [1] 0
## 
## $est.f.squared
## [1] 1.633993
# 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))
##      Comparison_i Comparison_j   mean_i   mean_j diff_in_means tukey.CI_lwr
## [1,]            1            2 29.12765 20.04512      9.082532     8.482181
## [2,]            1            3 29.12765 11.18160     17.946055    17.345704
## [3,]            2            3 20.04512 11.18160      8.863523     8.263171
##      tukey.CI_upr hedges_g        Q p        t           t_p  p.bonferroni
## [1,]     9.682883 1.586479 50.16889 0 35.47476 1.525250e-230 4.575749e-230
## [2,]    18.546406 3.134704 99.12803 0 70.09411  0.000000e+00  0.000000e+00
## [3,]     9.463874 1.548224 48.95915 0 34.61935 2.703667e-221 8.111001e-221
##      bonf.CI_lwr bonf.CI_upr
## [1,]    8.469263    9.695802
## [2,]   17.332785   18.559324
## [3,]    8.250253    9.476792
# t-test
lw.t.test(young, old, flr1 = min(young), ceil1 = max(young), flr2 = min(old), ceil2 = max(old), "a")
## $statistic
## [1] 69.32223
## 
## $p.value
## [1] 0
## 
## $est.d
## [1] 3.285657
## 
## $conf.int
## [1] 18.23623 19.29901
lw.t.test(young, old, flr1 = min(young), ceil1 = max(young), flr2 = min(old), ceil2 = max(old), "b")
## $statistic
## [1] 74.87169
## 
## $p.value
## [1] 0
## 
## $est.d
## [1] 3.348364
## 
## $conf.int
## [1] 18.27570 19.25954

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.

# 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)
## [1] 2.340214
## [1] 15.41184
mean(x.2);var(x.2)
## [1] 2.98688
## [1] 24.3313
# 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))
## $ceiling.percentage
## [1] 0
## 
## $floor.percentage
## [1] 0.23
## 
## $est.mean
## [1] 3.007903
## 
## $est.var
## [1] 4.604077
rec.mean.var(x.2.cf, flr = min(x.2.cf), ceil = max(x.2.cf))
## $ceiling.percentage
## [1] 0.09666667
## 
## $floor.percentage
## [1] 0
## 
## $est.mean
## [1] 2.978368
## 
## $est.var
## [1] 22.92892
# conduct a t test on healthy data
t.test(x.1,x.2)
## 
##  Welch Two Sample t-test
## 
## data:  x.1 and x.2
## t = -1.7767, df = 569.32, p-value = 0.07615
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.36156162  0.06823023
## sample estimates:
## mean of x mean of y 
##  2.340214  2.986880
t.test(x.1.cf,x.2.cf)
## 
##  Welch Two Sample t-test
## 
## data:  x.1.cf and x.2.cf
## t = 0.16229, df = 551.34, p-value = 0.8711
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.5788896  0.6831640
## sample estimates:
## mean of x mean of y 
##  2.780501  2.728364
# 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")
## $statistic
## [1] 0.09130877
## 
## $p.value
## [1] 0.9273421
## 
## $est.d
## [1] 0.007751838
## 
## $conf.int
## [1] -0.6084426  0.6675121
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")
## $statistic
## [1] 0.09749151
## 
## $p.value
## [1] 0.9224308
## 
## $est.d
## [1] 0.007960149
## 
## $conf.int
## [1] -0.5677351  0.6268045
# 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")
## $statistic
## [1] 1439.26
## 
## $p.value
## [1] 0
## 
## $est.f.squared
## [1] 0.09595066
lw.f.star(testdat.cf, y ~ group, flr = min(testdat.cf$y), ceil = max(testdat.cf$y), "b")
## $statistic
## [1] 1439.26
## 
## $p.value
## [1] 0
## 
## $est.f.squared
## [1] 0.09595066
# 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))
##      Comparison_i Comparison_j       mean_i     mean_j diff_in_means
## [1,]            1            2 -0.003270664  0.4067894    -0.4100601
## [2,]            1            3 -0.003270664 -0.3086038     0.3053331
## [3,]            2            3  0.406789388 -0.3086038     0.7153932
##      tukey.CI_lwr tukey.CI_upr   hedges_g        Q p         t           t_p
## [1,]   -0.4414187   -0.3787014 -0.4334183 43.34183 0 -30.64730 3.907763e-203
## [2,]    0.2739745    0.3366918  0.3227258 32.27258 0  22.82016 2.723277e-114
## [3,]    0.6840345    0.7467518  0.7561442 75.61442 0  53.46747  0.000000e+00
##       p.bonferroni bonf.CI_lwr bonf.CI_upr  gh.CI_lwr  gh.CI_upr games.howell.t
## [1,] 1.172329e-202  -0.4420932  -0.3780269 -0.4405961 -0.3795240      -31.47530
## [2,] 8.169830e-114   0.2732999   0.3373663  0.2721457  0.3385205       21.56429
## [3,]  0.000000e+00   0.6833600   0.7474264  0.6851050  0.7456814       55.36136
##      games.howell.p games.howell.g welch.satterthwaite.df welch.bonf.CI_lwr
## [1,]   1.358072e-08     -0.4333678               19246.30        -0.4412515
## [2,]   0.000000e+00      0.3226882               19993.23         0.2714334
## [3,]   1.065834e-08      0.7560560               19341.29         0.6844549
##      welch.bonf.CI_upr   t_Welch       p_Welch p_Welch.bonferroni
## [1,]        -0.3788686 -31.47530 4.381950e-212      1.314585e-211
## [2,]         0.3392328  21.56429 5.639017e-102      1.691705e-101
## [3,]         0.7463315  55.36136  0.000000e+00       0.000000e+00

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.

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")
## Step 1: Correcting marginal distributions...
##   Variables: 6  |  n: 400
##   Mean floor%: 7.0  |  Mean ceiling%: 6.3
## Step 2: Estimating pairwise correlations (MOM)...
##   Pair 3 / 15  Pair 6 / 15  Pair 9 / 15  Pair 12 / 15  Pair 15 / 15  Pair 15 / 15
## Step 3: nearPD correction: 0.00%
## Step 4: Estimating network...
## Note: Network with lowest lambda selected as best network: assumption of sparsity might be violated.
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret
## the presence of the smallest edges with care. Setting threshold = TRUE will
## enforce higher specificity, at the cost of sensitivity.
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret
## the presence of the smallest edges with care. Setting threshold = TRUE will
## enforce higher specificity, at the cost of sensitivity.
print(fit)
## 
## === GGM with Ceiling/Floor Correction (ggm_cfe) ===
## Variables: 6  |  Observations: 400
## Method(s): EBICglasso
## 
## Censoring summary (floor / ceiling %):
##   Mean:  7.0% / 6.3%
##   Range: 4.5-9.0% / 3.5-8.8%
## 
## nearPD correction: 0.00%
## 
## Estimated edges (corrected): 14  |  Density: 93.3%
## Estimated edges (naive):    14  |  Density: 93.3%
summary(fit)
## 
## === Censoring Statistics per Variable ===
##  variable pct_floor pct_ceil mu_naive mu_corrected sd_naive sd_corrected
##        V1       6.2      6.2 -0.01646     -0.01881    0.846        0.934
##        V2       8.2      3.5 -0.03577     -0.03719    0.838        0.912
##        V3       4.5      8.8  0.05091      0.05674    0.866        0.966
##        V4       7.2      4.2 -0.03208     -0.03960    0.892        1.012
##        V5       7.0      6.8 -0.02751     -0.03177    0.873        0.981
##        V6       9.0      8.5 -0.00148     -0.00132    0.909        1.048
## 
## nearPD correction: 0.000%
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.

Sigma.corrected=fit$Sigma_corrected

qgraph(Sigma.corrected,
       graph="glasso",
       sampleSize=fit$n,
       gamma=.5,
       layout="spring")
## Note: Network with lowest lambda selected as best network: assumption of sparsity might be violated.
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret
## the presence of the smallest edges with care. Setting threshold = TRUE will
## enforce higher specificity, at the cost of sensitivity.

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.

fit.all=ggm_cfe(y.cf, floor=-1.5, ceiling=1.5, method="all")
## Step 1: Correcting marginal distributions...
##   Variables: 6  |  n: 400
##   Mean floor%: 7.0  |  Mean ceiling%: 6.3
## Step 2: Estimating pairwise correlations (MOM)...
##   Pair 3 / 15  Pair 6 / 15  Pair 9 / 15  Pair 12 / 15  Pair 15 / 15  Pair 15 / 15
## Step 3: nearPD correction: 0.00%
## Step 4: Estimating network...
## Note: Network with lowest lambda selected as best network: assumption of sparsity might be violated.
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret
## the presence of the smallest edges with care. Setting threshold = TRUE will
## enforce higher specificity, at the cost of sensitivity.
## Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
## penalize.diagonal, : A dense regularized network was selected (lambda < 0.1 *
## lambda.max). Recent work indicates a possible drop in specificity. Interpret
## the presence of the smallest edges with care. Setting threshold = TRUE will
## enforce higher specificity, at the cost of sensitivity.
sapply(fit.all$networks_corrected, function(x) sum(x[upper.tri(x)]!=0))
## EBICglasso        FDR Bonferroni 
##         14         11          5
sapply(fit.all$networks_naive, function(x) sum(x[upper.tri(x)]!=0))
## EBICglasso        FDR Bonferroni 
##         14         11          4

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.