Consider a situation where we have multiple items of some sort which have associated effects, but we can’t association observations with single items; instead, each observation is associated with a group of items. (Two examples which have come up are (1) authorship of papers and (2) hockey players on a team.) These are not quite identical to “multi-membership models” (I think), although similar techniques to those shown here could work for multi-membership models.

Load packages (we don’t really need anything beyond lme4; the rest are for convenience/drawing pictures).

library(lme4)
library(broom.mixed)
library(ggplot2); theme_set(theme_bw())
library(Matrix)

Construct a simulated example: first, simulate the design (structure).

nm <- 20
nobs <- 500
set.seed(101)
## choose items for observations
pres <- matrix(rbinom(nobs*nm, prob=0.25, size=1),nrow=nobs,ncol=nm)
dimnames(pres) <- list(NULL,LETTERS[seq(nm)])
pres[1:5,]
##      A B C D E F G H I J K L M N O P Q R S T
## [1,] 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0
## [2,] 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## [3,] 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0
## [4,] 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0
## [5,] 0 1 0 1 0 0 0 0 0 0 1 0 1 1 1 1 1 0 0 0
table(rowSums(pres))
## 
##   1   2   3   4   5   6   7   8   9  10  11 
##  11  40  63  85 106  85  52  37  11   7   3

The first 10 observations:

image(Matrix(pres),ylim=c(1,10),sub="",ylab="Observation",
      xlab="Item",
      ## draw tick labels at top
      scales=list(at=1:20,x=list(labels=colnames(pres),
                                 alternating=2)))

Since we chose which items/individuals to include for each observation randomly and independently (Bernoulli with probability 0.25), we have a highly variable number of individuals present for different observations (1-11). This would be realistic for some examples (authorship), unrealistic for others (hockey) … I don’t think it really makes much difference computationally (statistically, having some observations with a single member must make estimation more powerful …) The 0/1 matrix (indicator variable for whether item \(i\) is included in observation \(j\)) is convenient, and will turn out to be the form we need for inclusion in the model. It should be fairly straightforward to convert other forms (e.g. a list of sets of items associated with each observation) to this form …

Now simulate the response variable.

b <- rnorm(nm)  ## item-level effects
## n.b. get in trouble if we don't add residual error
## (theta is scaled relative to residual error)
## here, use theta=sigma=1
y <- c(pres %*% b) +rnorm(nobs,sd=1)

Fit (see ?modular):

## helpful to specify a factor with the right levels:
## actual values are unimportant since we will specify Zt/Ztlist directly
fake <- rep(LETTERS[seq(nm)],length.out=nobs)
lmod <- lFormula(y~1+(1|fake))
lmod$reTrms$Zt <- lmod$reTrms$Ztlist[[1]] <- Matrix(t(pres))
devfun <- do.call(mkLmerDevfun, lmod)
opt <- optimizeLmer(devfun)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)

Results look OK (correct item and residual variance estimated):

m1
## Linear mixed model fit by REML ['lmerMod']
## REML criterion at convergence: 1500.701
## Random effects:
##  Groups   Name        Std.Dev.
##  fake     (Intercept) 0.9444  
##  Residual             0.9896  
## Number of obs: 500, groups:  fake, 20
## Fixed Effects:
## (Intercept)  
##      0.1392

The conditional modes by item look good:

dd <- tidy(m1, effects="ran_vals")
dd <- transform(dd, level=reorder(level,estimate))
truth <- data.frame(level=LETTERS[seq(nm)],estimate=b)
ggplot(dd,aes(x=level,y=estimate))+
    geom_pointrange(aes(ymin=estimate-2*std.error,
              ymax=estimate+2*std.error))+coord_flip()+
    geom_point(data=truth,colour="red")

Session info:

sessionInfo()
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] bindrcpp_0.2.2    ggplot2_3.1.0     broom.mixed_0.2.3 lme4_1.1-19      
## [5] Matrix_1.2-15    
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.0       plyr_1.8.4       glmmTMB_0.2.2.0  compiler_3.5.0  
##  [5] pillar_1.3.0     nloptr_1.2.1     bindr_0.1.1      TMB_1.7.15      
##  [9] tools_3.5.0      digest_0.6.18    gtable_0.2.0     evaluate_0.12   
## [13] tibble_1.4.2     nlme_3.1-137     lattice_0.20-38  pkgconfig_2.0.2 
## [17] rlang_0.3.0.1    yaml_2.2.0       coda_0.19-2      withr_2.1.2     
## [21] dplyr_0.7.8      stringr_1.3.1    knitr_1.20       rprojroot_1.3-2 
## [25] grid_3.5.0       tidyselect_0.2.5 glue_1.3.0       R6_2.3.0        
## [29] rmarkdown_1.10   minqa_1.2.4      reshape2_1.4.3   tidyr_0.8.2     
## [33] purrr_0.2.5      magrittr_1.5     scales_1.0.0     backports_1.1.2 
## [37] htmltools_0.3.6  MASS_7.3-51.1    splines_3.5.0    assertthat_0.2.0
## [41] colorspace_1.3-2 labeling_0.3     stringi_1.2.4    lazyeval_0.2.1  
## [45] munsell_0.5.0    broom_0.5.0      crayon_1.3.4