### Install lmerMultiMember from GitHub with devtools Source: https://jvparidon.github.io/lmerMultiMember/index.html Install the lmerMultiMember package from GitHub using devtools. The code first checks if devtools is installed and installs it if necessary, then installs the package from the repository jvparidon/lmerMultiMember. ```r if(!require(devtools)){ install.packages("devtools") library("devtools") } install_github("jvparidon/lmerMultiMember") ``` -------------------------------- ### Example: creating weight matrix from membership columns Source: https://jvparidon.github.io/lmerMultiMember/reference/weights_from_columns.html Creates a matrix a with membership columns and passes it to weights_from_columns to generate the weight matrix Wa. ```r a <- cbind( c("k", "k", "l", NA, "m"), c("l", "m", "k", "m", "k"), c("m", NA, "m", NA, NA) ) Wa <- weights_from_columns(a) ``` -------------------------------- ### Example: Create indicator matrix from sparse winner/loser matrices Source: https://jvparidon.github.io/lmerMultiMember/reference/bradleyterry_from_sparse.html Example usage creating sparse factor matrices from character vectors and passing them to bradleyterry_from_sparse. The result is stored in Wwl. ```r winners <- Matrix::fac2sparse(c("k", "k", "l", "m", "m")) losers <- Matrix::fac2sparse(c("l", "m", "m", "k", "l")) Wwl <- bradleyterry_from_sparse(winners, losers) ``` -------------------------------- ### Example: Creating indicator matrix from winner and loser vectors Source: https://jvparidon.github.io/lmerMultiMember/reference/bradleyterry_from_vectors.html Example usage of bradleyterry_from_vectors with sample winner and loser vectors. The resulting matrix Wwl is a sparse dgCMatrix. ```r winners <- c("k", "k", "l", "m", "m") losers <- c("l", "m", "m", "k", "l") Wwl <- bradleyterry_from_vectors(winners, losers) ``` -------------------------------- ### Example: creating interaction weights with interaction_weights() Source: https://jvparidon.github.io/lmerMultiMember/reference/interaction_weights.html Creates two sparse weight matrices using weights_from_vector() and Matrix::fac2sparse(), then computes their interaction with interaction_weights(). The result Wab is a sparse interaction weight/indicator matrix of class Matrix::dgCMatrix. ```r a <- rep(c("k", "l", "l,k"), 2) b <- rep(c("m", "n"), 3) Wa <- weights_from_vector(a) Wb <- Matrix::fac2sparse(b) Wab <- interaction_weights(Wa, Wb) ``` -------------------------------- ### Example: fitting a multiple membership model Source: https://jvparidon.github.io/lmerMultiMember/reference/lmer.html Creates a data frame with a memberships column, converts it to weights, and fits a model with a custom membership matrix. The output shows the model fit and estimates. ```r df <- data.frame( x = seq(60) + runif(60, 0, 10), y = seq(60) + rep(runif(6, 0, 10), 10), memberships = rep(c("a,b,c", "a,c", "a", "b", "b,a", "b,c,a"), 10) ) weights <- weights_from_vector(df$memberships) # note that the grouping variable name is arbitrary -- it just has # to match the name in the list and doesn't need to correspond to a column # name in the data lmer(y ~ x + (1 | members), data = df, memberships = list(members = weights) ) #> Linear mixed model fit by REML ['lmerModMultiMember'] #> Formula: y ~ x + (1 | members) #> Data: df #> REML criterion at convergence: 334.1473 #> Random effects: #> Groups Name Std.Dev. #> members (Intercept) 4.081e-07 #> Residual 3.827e+00 #> Number of obs: 60, groups: members, 3 #> Fixed Effects: #> (Intercept) x #> -0.2640 0.9876 ``` -------------------------------- ### Load packages for group membership example Source: https://jvparidon.github.io/lmerMultiMember/articles/group_membership.html Loads lmerMultiMember and convenience packages (broom.mixed, ggplot2, Matrix) for modeling and plotting. Only lmerMultiMember is strictly required. ```r library(lmerMultiMember) library(broom.mixed) library(ggplot2) theme_set(theme_bw()) library(Matrix) ``` -------------------------------- ### weights_from_vector example Source: https://jvparidon.github.io/lmerMultiMember/reference/weights_from_vector.html Creates a membership vector a with comma-separated group memberships and calls weights_from_vector(a) to produce the sparse weight matrix Wa. The result is a Matrix::dgCMatrix object. ```r a <- c("k,l,m", "k,m", "k", "l", "l,k") Wa <- weights_from_vector(a) ``` -------------------------------- ### Load R packages and display NFL data Source: https://jvparidon.github.io/lmerMultiMember/articles/bradleyterry_models.html Initializes the environment with required libraries and renders the nfl_scores_2021 dataset as a scrollable, styled HTML table. ```R # load packages library(lmerMultiMember) # for dataset and modeling library(kableExtra) # for displaying tabular data library(dplyr) # for data manipulation library(sjPlot) # for plotting model coefficients library(ggplot2) # for general plotting stuff # display the NFL 2021 season data nfl_scores_2021 |> kable() |> kable_styling(bootstrap_options = c("hover", "condensed"), font_size = 16) |> scroll_box(height = "250px") ``` -------------------------------- ### Set up Stan sampling and prepare data Source: https://jvparidon.github.io/lmerMultiMember/articles/bradleyterry_models.html Configures multiple cores for sampling and constructs the data list required by the Stan model, including the home win indicator, number of observations, number of teams, and the design matrix from the bradleyterry_from_vectors function. ```R # use multiple cores for sampling options(mc.cores = parallel::detectCores()) # prepare data for passing to Stan sampler stan_data <- list( home_win = nfl_scores_2021$home_win, N = nrow(nfl_scores_2021), M = length(unique(nfl_scores_2021$home_team)), W = t(as.matrix(bradleyterry_from_vectors(nfl_scores_2021$home_team, nfl_scores_2021$visiting_team))) ) # sample from Stan model stan_fit <- rstan::sampling(model_stan, data = stan_data, warmup = 2000, iter = 8000) # create tidy data frame from posterior samples stan_post <- broom.mixed::tidy(stan_fit, pars = "team_OR", conf.int = TRUE) # store fitted Stan model save(stan_post, file = "bt_model_stan_fit.Rds") ``` -------------------------------- ### glmer(formula, data, family, control, start, verbose, nAGQ, weights, na.action, offset, contrasts, devFunOnly, memberships) Source: https://jvparidon.github.io/lmerMultiMember/reference/glmer.html Fits a generalized linear mixed effects model with multimembership random effects. This function extends lme4::glmer by accepting a memberships argument to specify multimembership structures. ```APIDOC ## glmer(formula, data, family, control, start, verbose, nAGQ, weights, na.action, offset, contrasts, devFunOnly, memberships) ### Description Fits a generalized linear mixed effects model with multimembership random effects. This function extends lme4::glmer by accepting a memberships argument to specify multimembership structures. ### Usage ```r glmer( formula, data = NULL, family, control = lme4::glmerControl(), start = NULL, verbose = 0L, nAGQ = 1L, weights = NULL, na.action = na.omit, offset = NULL, contrasts = NULL, devFunOnly = FALSE, memberships = NULL ) ``` ### Arguments - **formula** (formula) - Required - A two-sided linear formula object describing the fixed and random effects parts of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right. - **data** (data.frame) - Optional - An optional data frame containing the variables named in formula. By default the variables are taken from the environment from which glmer is called. - **family** (family) - Required - A GLM family, see glm and family for details. - **control** (list) - Optional - A list of control options, see lme4::glmerControl for details. - **start** (list) - Optional - A list of starting values for the parameters, see lme4::glmer for details. - **verbose** (integer) - Optional - Integer scalar indicating the level of output during fitting. Default is 0L. - **nAGQ** (integer) - Optional - Integer scalar - the number of points per axis for evaluating the adaptive Gauss-Hermite approximation to the log-likelihood. Default is 1L. - **weights** (vector) - Optional - An optional vector of prior weights to be used in the fitting process. Should be NULL or a numeric vector. - **na.action** (function) - Optional - A function that indicates what should happen when the data contain NA values. Default is na.omit. - **offset** (numeric) - Optional - A numeric vector of length equal to the number of cases to be used as an offset in the model. - **contrasts** (list) - Optional - An optional list of contrasts to be used for the fixed effect terms. - **devFunOnly** (logical) - Optional - Logical indicating whether to return the deviance function only. Default is FALSE. - **memberships** (list) - Optional - A list specifying the multimembership structure for the random effects. See package documentation for details. ### Details This function is a wrapper around lme4::glmer that allows for multimembership random effects. The memberships argument is used to define how observations belong to multiple random effect levels. ### Value An object of class merMod, similar to the output of lme4::glmer, but with multimembership random effects incorporated. ### References lme4::glmer documentation. ``` -------------------------------- ### lmer(formula, data, REML, control, start, verbose, weights, na.action, offset, contrasts, devFunOnly, memberships) Source: https://jvparidon.github.io/lmerMultiMember/reference/lmer.html Fits a linear mixed-effects model with optional multi-member random effects. The memberships argument accepts a named list of weight matrices that replace dummy random effects with matching names. ```APIDOC ## lmer(formula, data, REML, control, start, verbose, weights, na.action, offset, contrasts, devFunOnly, memberships) ### Description Fits a linear mixed-effects model with optional multi-member random effects. The memberships argument accepts a named list of weight matrices that replace dummy random effects with matching names. ### Usage ```r lmer( formula, data = NULL, REML = TRUE, control = lme4::lmerControl(), start = NULL, verbose = 0L, weights = NULL, na.action = na.omit, offset = NULL, contrasts = NULL, devFunOnly = FALSE, memberships = NULL ) ``` ### Arguments - **formula** (formula) - Required - A two-sided linear formula object describing both the fixed-effects and random-effects part of the model, with the response on the left of a `~` operator and the terms, separated by `+` operators, on the right. Random-effects terms are distinguished by vertical bars (`|`) separating expressions for design matrices from grouping factors. Two vertical bars (`||`) can be used to specify multiple uncorrelated random effects for the same grouping variable. (Because of the way it is implemented, the `||`-syntax _works only for design matrices containing numeric (continuous) predictors_ ; to fit models with independent categorical effects, see `dummy` or the `lmer_alt` function from the afex package.) - **data** (data.frame) - Optional - An optional data frame containing the variables named in `formula`. By default the variables are taken from the environment from which `lmer` is called. While `data` is optional, the package authors _strongly_ recommend its use, especially when later applying methods such as `update` and `drop1` to the fitted model (_such methods are not guaranteed to work properly if`data` is omitted_). If `data` is omitted, variables will be taken from the environment of `formula` (if specified as a formula) or from the parent frame (if specified as a character vector). - **REML** (logical) - Optional - Logical scalar - Should the estimates be chosen to optimize the REML criterion (as opposed to the log-likelihood)? - **control** (list) - Optional - A list (of correct class, resulting from `lmerControl()` or `glmerControl()` respectively) containing control parameters, including the nonlinear optimizer to be used and parameters to be passed through to the nonlinear optimizer, see the `*lmerControl` documentation for details. - **start** (list) - Optional - A named `list` of starting values for the parameters in the model. For `lmer` this can be a numeric vector or a list with one component named `"theta"`. - **verbose** (integer) - Optional - Integer scalar. If `> 0` verbose output is generated during the optimization of the parameter estimates. If `> 1` verbose output is generated during the individual penalized iteratively reweighted least squares (PIRLS) steps. - **weights** (numeric vector) - Optional - An optional vector of ‘prior weights’ to be used in the fitting process. Should be `NULL` or a numeric vector. Prior `weights` are _not_ normalized or standardized in any way. In particular, the diagonal of the residual covariance matrix is the squared residual standard deviation parameter `sigma` times the vector of inverse `weights`. Therefore, if the `weights` have relatively large magnitudes, then in order to compensate, the `sigma` parameter will also need to have a relatively large magnitude. - **na.action** (function) - Optional - A function that indicates what should happen when the data contain `NA`s. The default action (`na.omit`, inherited from the 'factory fresh' value of `getOption("na.action")`) strips any observations with any missing values in any variables. - **offset** (numeric vector) - Optional - This can be used to specify an _a priori_ known component to be included in the linear predictor during fitting. This should be `NULL` or a numeric vector of length equal to the number of cases. One or more `offset` terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See `model.offset`. - **contrasts** (list) - Optional - An optional list. See the `contrasts.arg` of `model.matrix.default`. - **devFunOnly** (logical) - Optional - Logical - return only the deviance evaluation function. Note that because the deviance function operates on variables stored in its environment, it may not return _exactly_ the same values on subsequent calls (but the results should always be within machine tolerance). - **memberships** (named list) - Optional - Named list of weight matrices that will replace any (dummy) random effects with matching names. ### Value lme4 model object ### Example ```r df <- data.frame( x = seq(60) + runif(60, 0, 10), y = seq(60) + rep(runif(6, 0, 10), 10), memberships = rep(c("a,b,c", "a,c", "a", "b", "b,a", "b,c,a"), 10) ) weights <- weights_from_vector(df$memberships) # note that the grouping variable name is arbitrary -- it just has # to match the name in the list and doesn't need to correspond to a column # name in the data lmer(y ~ x + (1 | members), data = df, memberships = list(members = weights) ) #> Linear mixed model fit by REML ['lmerModMultiMember'] #> Formula: y ~ x + (1 | members) #> Data: df #> REML criterion at convergence: 334.1473 #> Random effects: #> Groups Name Std.Dev. #> members (Intercept) 4.081e-07 #> Residual 3.827e+00 #> Number of obs: 60, groups: members, 3 #> Fixed Effects: #> (Intercept) x #> -0.2640 0.9876 ``` ``` -------------------------------- ### Visualize presence matrix with image() Source: https://jvparidon.github.io/lmerMultiMember/articles/group_membership.html Plots the simulated presence matrix using image() from the Matrix package, showing the first 10 observations and labeling axes as Observation and Item. ```r image(Matrix(pres), ylim = c(1, 10), sub = "", ylab = "Observation", xlab = "Item", ``` -------------------------------- ### Load ATP doubles data and display first 1000 matches Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Loads the atp_doubles dataset from the lmerMultiMember package and displays the first 1000 matches in a scrollable table using kableExtra. ```R # load ATP doubles data from package data("atp_doubles", package = "lmerMultiMember") # display the first 1000 matches in the ATP doubles dataset atp_doubles %>% head(1000) %>% kable() %>% kable_styling(bootstrap_options = c("hover", "condensed"), font_size = 16) %>% scroll_box(height = "250px") ``` -------------------------------- ### Generate and plot sparse indicator matrix for small dataset Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Creates a sparse indicator matrix for the first five matches using the weights_from_columns function, then plots it with plot_membership_matrix and updates the plot labels and colors. The matrix shows 1s for players in team 1 and 0s elsewhere. ```R # generate sparse matrix for small bit of the dataset Wt <- weights_from_columns(atp_doubles[1:5, c("team1_player1_name", "team1_player2_name")]) # display matrix plot of the indicator matrix matplot <- plot_membership_matrix(Wt) update(matplot, xlab = "Team 1 players", ylab = "Match", at = c(-1.001, .5, 1.001), col.regions = c("white", "coral1")) ``` -------------------------------- ### Generate Bradley-Terry indicator matrix and plot membership matrix Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Creates indicator matrices for team 1 and team 2 players from the first five rows of the atp_doubles dataset, combines them into a Bradley-Terry/adversarial indicator matrix using bradleyterry_from_sparse, and plots the membership matrix with a custom color palette. The plot uses blue for negative indicators and coral for positive ones. ```r # generate indicator matrices for both teams for part of the dataset Wp1 <- weights_from_columns(atp_doubles[1:5, c("team1_player1_name", "team1_player2_name")]) Wp2 <- weights_from_columns(atp_doubles[1:5, c("team2_player1_name", "team2_player2_name")]) # create a Bradley-Terry/adversarial indicator matrix Wp <- bradleyterry_from_sparse(Wp1, Wp2) # display matrix plot of the indicator matrix matplot <- plot_membership_matrix(Wp) cmap <- colorRampPalette(c("#6495ed", "white", "coral1"))(100) matplot$legend$right$args$key$col <- cmap update(matplot, xlab = "Players", ylab = "Match", at = seq(-1.001, 1.001, length.out = 100), col.regions = cmap) ``` -------------------------------- ### weights_from_columns usage Source: https://jvparidon.github.io/lmerMultiMember/reference/weights_from_columns.html Shows the function signature for weights_from_columns, which takes membership_columns and returns a sparse weight/indicator matrix of type Matrix::dgCMatrix. ```r weights_from_columns(membership_columns) ``` -------------------------------- ### Fit multiple membership GLMM and display summary Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Generates the full sparse indicator matrix for all matches, fits a binomial GLMM with a multiple membership random effect for team1_player, and displays the model summary. The summary shows the model is a multiple membership model with 2 players per match. ```R # generate sparse indicator matrix for whole dataset Wp <- weights_from_columns(atp_doubles[, c("team1_player1_name", "team1_player2_name")]) # fit GLMM m1 <- glmer(team1_win ~ 0 + (1 | team1_player), family = "binomial", memberships = list(team1_player = Wp), data = atp_doubles) # display model summary summary(m1) ``` ```R ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation). Model includes multiple membership random effects. [ ## glmerModMultiMember] ## Family: binomial ( logit ) ## Formula: team1_win ~ 0 + (1 | team1_player) ## Data: atp_doubles ## ## AIC BIC logLik deviance df.resid ## 17742.3 17749.8 -8870.1 17740.3 13037 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -1.71911 -0.99037 0.04637 0.91606 1.75388 ## ## Random effects: ## Groups Name Variance Std.Dev. ## team1_player (Intercept) 0.09373 0.3062 ## Number of obs: 13038, groups: team1_player, 858 ## ## Group memberships per observation for multiple membership REs: ## Min. per obs. Mean per obs. Max. per obs. ## team1_player 2 2 2 ``` -------------------------------- ### Fit multiple membership model with player-year interaction Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Creates a sparse interaction matrix from player and year factors, fits a binomial glmer with multiple membership random effects, and prints the model summary. The model uses the formula team1_win ~ 0 + (1 | playerXyear) with the memberships argument specifying the interaction matrix. ```R # get year in which tournament was played atp_doubles$year <- substr(atp_doubles$tourney_date, 1, 4) # create sparse matrix from years Wy <- Matrix::fac2sparse(atp_doubles$year) # generate sparse interaction matrix from sparse player and year matrices Wpy <- interaction_weights(Wp, Wy) # fit model m3 <- glmer(team1_win ~ 0 + (1 | playerXyear), family = "binomial", memberships = list(playerXyear = Wpy), data = atp_doubles) # print model summary summary(m3) ``` ```R ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation). Model includes multiple membership random effects. [ ## glmerModMultiMember] ## Family: binomial ( logit ) ## Formula: team1_win ~ 0 + (1 | playerXyear) ## Data: atp_doubles ## ## AIC BIC logLik deviance df.resid ## 17297.1 17304.6 -8647.5 17295.1 13037 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -3.14747 -0.84233 -0.01485 0.85261 2.37473 ## ## Random effects: ## Groups Name Variance Std.Dev. ## playerXyear (Intercept) 0.2066 0.4546 ## Number of obs: 13038, groups: playerXyear, 9920 ## ## Group memberships per observation for multiple membership REs: ## Min. per obs. Mean per obs. Max. per obs. ## playerXyear 0 0 0 ``` -------------------------------- ### Display model summary with summary(m2) Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Displays the summary of the Bradley-Terry style multiple membership model m2. The output shows the model formula, fit statistics, random effects variance, and a multiple membership summary line. The summary indicates that the indicators for team 1 and team 2 cancel out, resulting in min/mean/max of 0 for observations, which is expected for a correctly specified Bradley-Terry model. ```r summary(m2) ``` ```output ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation). Model includes multiple membership random effects. [ ## glmerModMultiMember] ## Family: binomial ( logit ) ## Formula: team1_win ~ 0 + (1 | player) ## Data: atp_doubles ## ## AIC BIC logLik deviance df.resid ## 16688.3 16695.8 -8343.1 16686.3 13037 ## ## Scaled residuals: ## Min 1Q Median 3Q Max ## -2.98010 -0.84939 0.00367 0.85098 2.85597 ## ## Random effects: ## Groups Name Variance Std.Dev. ## player (Intercept) 0.1903 0.4363 ## Number of obs: 13038, groups: player, 992 ## ## Group memberships per observation for multiple membership REs: ## Min. per obs. Mean per obs. Max. per obs. ## player 0 0 0 ``` -------------------------------- ### Session information Source: https://jvparidon.github.io/lmerMultiMember/articles/group_membership.html Output of sessionInfo() showing the R version, platform, locale, and attached packages used in the analysis. ```R ## R version 4.3.1 (2023-06-16) ## Platform: x86_64-pc-linux-gnu (64-bit) ## Running under: Ubuntu 22.04.3 LTS ## ## Matrix products: default ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0 ## ## locale: ## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 ## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 ## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C ## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C ## ## time zone: UTC ## tzcode source: system (glibc) ## ## attached base packages: ## [1] stats graphics grDevices utils datasets methods base ## ## other attached packages: ## [1] Matrix_1.5-4.1 ggplot2_3.4.4 broom.mixed_0.2.9.4 ## [4] lmerMultiMember_0.11.8 ## ## loaded via a namespace (and not attached): ## [1] future_1.33.0 sass_0.4.7 utf8_1.2.3 generics_0.1.3 ## [5] tidyr_1.3.0 stringi_1.7.12 lattice_0.21-8 listenv_0.9.0 ## [9] lme4_1.1-34 digest_0.6.33 magrittr_2.0.3 evaluate_0.22 ## [13] grid_4.3.1 fastmap_1.1.1 rprojroot_2.0.3 jsonlite_1.8.7 ## [17] backports_1.4.1 purrr_1.0.2 fansi_1.0.5 scales_1.2.1 ## [21] codetools_0.2-19 textshaping_0.3.7 jquerylib_0.1.4 cli_3.6.1 ## [25] rlang_1.1.1 parallelly_1.36.0 munsell_0.5.0 splines_4.3.1 ## [29] withr_2.5.1 cachem_1.0.8 yaml_2.3.7 parallel_4.3.1 ## [33] tools_4.3.1 memoise_2.0.1 nloptr_2.0.3 minqa_1.2.6 ## [37] dplyr_1.1.3 colorspace_2.1-0 forcats_1.0.0 globals_0.16.2 ## [41] boot_1.3-28.1 broom_1.0.5 vctrs_0.6.4 R6_2.5.1 ## [45] lifecycle_1.0.3 stringr_1.5.0 fs_1.6.3 MASS_7.3-60 ## [49] furrr_0.3.1 ragg_1.2.6 pkgconfig_2.0.3 desc_1.4.2 ## [53] gtable_0.3.4 pkgdown_2.0.7 bslib_0.5.1 pillar_1.9.0 ## [57] glue_1.6.2 Rcpp_1.0.11 systemfonts_1.0.5 xfun_0.40 ## [61] tibble_3.2.1 tidyselect_1.2.0 knitr_1.44 farver_2.1.1 ## [65] htmltools_0.5.6.1 nlme_3.1-162 labeling_0.4.3 rmarkdown_2.25 ## [69] compiler_4.3.1 ``` -------------------------------- ### Simulate binary presence matrix for group membership Source: https://jvparidon.github.io/lmerMultiMember/articles/group_membership.html Creates a simulated design with 20 items and 1000 observations, where each observation has a 25% chance of including each item. The first five rows are printed to show the binary presence matrix. ```r nm <- 20 # number of groups nobs <- 1000 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, ] ``` ```text ## 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 1 0 0 0 0 0 1 0 0 1 1 1 1 0 0 0 ## [2,] 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 1 ## [3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 ## [4,] 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 1 0 0 ## [5,] 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 ``` -------------------------------- ### Usage: bradleyterry_from_sparse(winners, losers) Source: https://jvparidon.github.io/lmerMultiMember/reference/bradleyterry_from_sparse.html Function signature for bradleyterry_from_sparse. Takes two sparse matrices: winners and losers. Returns a sparse Bradley-Terry indicator matrix of type Matrix::dgCMatrix. ```r bradleyterry_from_sparse(winners, losers) ``` -------------------------------- ### interaction_weights usage Source: https://jvparidon.github.io/lmerMultiMember/reference/interaction_weights.html Shows the function signature for interaction_weights(). It accepts two sparse weight/indicator matrices of class Matrix::dgCMatrix. ```r interaction_weights(a, b) ``` -------------------------------- ### Generate Bradley-Terry indicator matrix with bradleyterry_from_vectors() Source: https://jvparidon.github.io/lmerMultiMember/articles/bradleyterry_models.html Use this snippet to create a Bradley-Terry indicator matrix from the home and visiting team vectors in the NFL 2021 dataset. The function is from the lmerMultiMember package and takes the home_team and visiting_team columns as arguments. ```R W <- bradleyterry_from_vectors(nfl_scores_2021$home_team, nfl_scores_2021$visiting_team) ``` -------------------------------- ### Simulated data head output Source: https://jvparidon.github.io/lmerMultiMember/articles/group_membership.html First six rows of the simulated dataset showing columns x, z, g1, g2, grps, and y. ```r ## x z g1 g2 grps y ## 1 0.9148060 0.5283896 11 13 K,M 9.7194660 ## 2 0.9370754 0.6463879 24 12 X,L 4.9175119 ## 3 0.2861395 0.8340490 17 2 Q,B 0.3754283 ## 4 0.8304476 0.3457626 22 19 V,S 8.4488474 ## 5 0.6417455 0.6217329 25 20 Y,T 5.8487113 ## 6 0.5190959 0.6631953 23 10 W,J 9.5083034 ``` -------------------------------- ### weights_from_columns(membership_columns) Source: https://jvparidon.github.io/lmerMultiMember/reference/weights_from_columns.html Generates a sparse weight/indicator matrix from columns containing group memberships. The input should be a matrix or data frame where each column represents a group membership variable, and the output is a sparse matrix of class dgCMatrix. ```APIDOC ## weights_from_columns(membership_columns) ### Description Generates a sparse weight/indicator matrix from columns containing group memberships. ### Usage ```r weights_from_columns(membership_columns) ``` ### Arguments - **membership_columns** - columns containing group memberships ### Value sparse weight/indicator matrix of type Matrix::dgCMatrix ### Examples ```r a <- cbind( c("k", "k", "l", NA, "m"), c("l", "m", "k", "m", "k"), c("m", NA, "m", NA, NA) ) Wa <- weights_from_columns(a) ``` ``` -------------------------------- ### print(x, ...) Source: https://jvparidon.github.io/lmerMultiMember/reference/print.summary.lmerModLmerTestMultiMember.html S3 method for summary.lmerModLmerTestMultiMember objects. Prints the summary of a multiple memberships model, optionally passing additional arguments to print.summary.merMod. ```APIDOC ## print(x, ...) ### Description S3 method for printing a summary.lmerModLmerTestMultiMember object, which is the summary of a multiple memberships model fitted with lmerMultiMember. ### Usage ```r # S3 method for summary.lmerModLmerTestMultiMember print(x, ...) ``` ### Arguments - **x** (merModMultiMember object) - Required - The summary object to print. - **...** - Optional - additional arguments to be passed on to print.summary.merMod. ### Details This method allows users to display the summary of a multiple memberships model in a formatted way. The additional arguments are forwarded to the standard print method for merMod summaries. ``` -------------------------------- ### Simulate data with multiple membership and random slopes Source: https://jvparidon.github.io/lmerMultiMember/articles/group_membership.html Simulates a dataset with two grouping factors (g1, g2) and random slopes for x, including fixed effects and random effects. The group contributions are assumed independent and equally weighted. ```r set.seed(42) nobs <- 10000 nm <- 26 dat <- data.frame( x = runif(nobs), z = runif(nobs), # TODO: double check that all(g1 != g2) g1 = sample(1:nm, nobs, TRUE), g2 = sample(1:nm, nobs, TRUE) ) dat$grps <- paste(LETTERS[dat$g1], LETTERS[dat$g2], sep = ",") # fixed effects beta <- c(1, 2, 3) # random effects bint <- rnorm(nm, sd = 3.14) bx <- rnorm(nm, sd = 0.2) # we could be clever about setting the model matrix here, but it's # easier and clearer how we're constructing bits if we're inefficient fe <- beta[1] + beta[2] * dat$x + beta[3] * dat$z # note that the group contributions for each observation are assumed to be # independent and equally weighted. # TODO: add example of unequal weights re <- bint[dat$g1] + bint[dat$g2] + (bx[dat$g1] + bx[dat$g2]) * dat$x dat$y <- fe + re + rnorm(nobs, sd = 1) head(dat) ``` -------------------------------- ### Fit multi-membership GLMM with glmer() and memberships argument Source: https://jvparidon.github.io/lmerMultiMember/reference/glmer.html Fits a multi-membership GLMM with glmer() by passing a named list of weight matrices to the memberships argument. The grouping variable name in the formula is arbitrary and only needs to match the list name, not a column in the data. The displayed output shows the resulting glmerModMultiMember summary. ```r df <- data.frame( x = runif(60, 0, 1), y = rbinom(60, 1, 0.6), memberships = rep(c("a,b,c", "a,c", "a", "b", "b,a", "b,c,a"), 10) ) weights <- weights_from_vector(df$memberships) # note that the grouping variable name is arbitrary -- it just has # to match the name in the list and doesn't need to correspond to a column # name in the data glmer(y ~ x + (1 | members), data = df, family = binomial, memberships = list(members = weights) ) ``` ```r #> Generalized linear mixed model fit by maximum likelihood (Laplace #> Approximation) [glmerModMultiMember] #> Family: binomial ( logit ) #> Formula: y ~ x + (1 | members) #> Data: df #> AIC BIC logLik deviance df.resid #> 83.4162 89.6992 -38.7081 77.4162 57 #> Random effects: #> Groups Name Std.Dev. #> members (Intercept) 0 #> Number of obs: 60, groups: members, 3 #> Fixed Effects: #> (Intercept) x #> 0.3690 0.4972 ``` -------------------------------- ### interaction_weights(a, b) Source: https://jvparidon.github.io/lmerMultiMember/reference/interaction_weights.html Creates a sparse interaction weight/indicator matrix from two input sparse matrices of class dgCMatrix. ```APIDOC ## interaction_weights(a, b) ### Description Provides a helper function for pre-generating interactions. Takes sparse weight/indicator matrices generated with e.g. Matrix::fac2sparse() or weights_from_vector() as input. ### Usage ```r interaction_weights(a, b) ``` ### Arguments - **a** (Matrix::dgCMatrix) - Required - sparse weight/indicator matrix of class Matrix::dgCMatrix - **b** (Matrix::dgCMatrix) - Required - sparse weight/indicator matrix of class Matrix::dgCMatrix ### Value sparse interaction weight/indicator matrix of class Matrix::dgCMatrix ### Examples ```r a <- rep(c("k", "l", "l,k"), 2) b <- rep(c("m", "n"), 3) Wa <- weights_from_vector(a) Wb <- Matrix::fac2sparse(b) Wab <- interaction_weights(Wa, Wb) ``` ``` -------------------------------- ### ANOVA comparison with anova(m1, m2) Source: https://jvparidon.github.io/lmerMultiMember/articles/lmermultimember_intro.html Performs an ANOVA comparison of models m1 and m2. The output shows that model 2 has a lower log-likelihood and deviance with the same number of parameters (Df difference = 0), confirming that the Bradley-Terry style model fits better without adding parameters. The p-value is highly significant. ```r anova(m1, m2) ``` ```output ## Data: atp_doubles ## Models: ## m1: team1_win ~ 0 + (1 | team1_player) ## m2: team1_win ~ 0 + (1 | player) ## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq) ## m1 1 17742 17750 -8870.1 17740 ## m2 1 16688 16696 -8343.1 16686 1054 0 ``` -------------------------------- ### plot_membership_matrix.summary.merModMultiMember Source: https://jvparidon.github.io/lmerMultiMember/reference/plot_membership_matrix.summary.merModMultiMember.html Documentation for the S3 method plot_membership_matrix for summary.merModMultiMember objects. ```APIDOC ## plot_membership_matrix.summary.merModMultiMember ### Description Plots the membership matrix for a multiple memberships model summary object. ### Usage ```r # S3 method for summary.merModMultiMember plot_membership_matrix(x, varname, ...) ``` ### Arguments - **x** (merModMultiMember summary object) - Required - A summary object from a multiple memberships model. - **varname** (character) - Required - The multimembership variable to be plotted. - **...** - Optional - Additional arguments to pass on to `Matrix::image`. ### Details This is an S3 method for the `summary.merModMultiMember` class. It plots the membership matrix for the specified multimembership variable. ### Value No return value; used for its side effect of producing a plot. ``` -------------------------------- ### plot_membership_matrix(x, varname = NULL, ...) Source: https://jvparidon.github.io/lmerMultiMember/reference/plot_membership_matrix.default.html Default S3 method for plotting a multiple membership weight matrix. It creates a plot of the matrix, using the provided variable name as the title, and forwards additional arguments to Matrix::image. ```APIDOC ## plot_membership_matrix(x, varname = NULL, ...) ### Description Default S3 method for plotting a multiple membership weight matrix. It creates a plot of the matrix, using the provided variable name as the title, and forwards additional arguments to Matrix::image. ### Usage ```r # S3 method for default plot_membership_matrix(x, varname = NULL, ...) ``` ### Arguments - **x** - multimembership weight matrix - **varname** - name of the multimembership variable, will be used as title - **...** - additional arguments to pass on to Matrix::image ``` -------------------------------- ### Usage: bradleyterry_from_vectors(winners, losers) Source: https://jvparidon.github.io/lmerMultiMember/reference/bradleyterry_from_vectors.html Function signature for bradleyterry_from_vectors. It takes two vectors of strings: winners and losers, and returns a sparse Bradley-Terry indicator matrix of type Matrix::dgCMatrix. ```r bradleyterry_from_vectors(winners, losers) ``` -------------------------------- ### weights_from_vector usage Source: https://jvparidon.github.io/lmerMultiMember/reference/weights_from_vector.html Shows the function signature with the membership_vector and sep arguments. The sep argument defaults to "," and delimits group memberships in the input vector. ```r weights_from_vector(membership_vector, sep = ",") ```