Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
113,080 characters · 17 sections · 112 citation commands
Recidivism and Peer Influence with LLM Text Embeddings in Low Security Correctional Facilities
\spacingset{1.5}
JEL Classification: C31, C36, C55, K42.\\ Keywords: Peer effects, Language embeddings, Network endogeneity, Recidivism.
Peer effects play an important role in many outcomes economists are interested in, including educational attainment, criminal behavior, and workplace productivity bayer2009building,mas2009peers,sacerdote2001peer. Researchers have identified several mechanisms by which peers influence individuals' behavior boucher2024toward,bursztyn2014understanding. These include social learning, where individuals learn from or update beliefs based on peers' actions banerjee1992simple,bikhchandani1992theory, and models of social image and peer pressure, where individuals modify behavior to adapt to societal norms bursztyn2017social. In this paper, we focus on language as an underlying mechanism and propose studying peer effects on language. It is well known that individuals often convey behavioral traits and beliefs through language exchanges fouka2020backlash,grogger2019speech,giles2013communication,niederhoffer2002linguistic. Therefore, peers' language reveals their social alignment and beliefs, making it central to these mechanisms of peer effects.
There are two major challenges to estimating peer effects in language. First, language interactions, by nature, are non-numerical, unstructured, high-dimensional, and context-dependent, making them difficult to analyze mathematically. Second, similar to peer effect estimation in other settings, the identification of peer effects in language is complex due to multiple issues. These include simultaneity in the outcomes and confounding due to unobserved homophily affecting both the outcome and peer selection manski1993identification,shalizi2011homophily,johnsson2021estimation.
The first challenge is to represent language interactions in a structured numerical form. Traditional text analysis methods, such as bag-of-words, dictionary methods, and topic models have limitations in capturing implicit semantic meaning and often require explicit appearance of keywords specified by the researcher in the text. Earlier text embedding models, such as Word2vec mikolov2013efficient attempt to address these limitations by representing text as high-dimensional embedding vectors. However, even these models are limited by their inability to incorporate context in their embeddings. In terms of the second challenge, identification of peer effects can be achieved using instrumental variables, where peers interact in network settings, using the identification results outlined in bramoulle2009identification. However, bramoulle2009identification's results assume that the network is uncorrelated with the error term conditional on covariates. Recently, johnsson2021estimation have extended the results of bramoulle2009identification, allowing for network endogeneity. However, these results rely on strong assumptions about network density that are unrealistic for real-world network data. Consequently, the identification and estimation of peer effects in language require new econometric methods and theory that accommodate multivariate, correlated outcomes, endogeneity in network formation, and are developed under the assumption of realistic, sparse networks.
In this paper, we provide a unified framework that combines Large Language Model embeddings with novel econometric methodology for the identification and estimation of peer effects in language interactions. We use embeddings from transformer-based foundation AI models that improve upon early text embedding models by incorporating both an attention mechanism vaswani2017attention and bidirectional processing devlin2019bert to obtain context-aware vector representations of text. These models can generate different vector embeddings for the same word because other words in the sentence dynamically influence its embeddings, thereby defining a context for the word. Foundation models trained on extensive text corpora have demonstrated strong capabilities for learning language representations. Consequently, these models can be applied to obtain embeddings for new text in specific applications without the need for any new training Radford2018ImprovingLU,brown2020language.
Even though the LLM embedding vectors provide a meaningful numerical representation of language, they are not well-suited for peer effect estimation due to their high dimensionality (700-1000). Such high dimensionality limits the scope of interpretation in peer effect estimation. We address this problem using zero-shot classification, which provides a suitable way to map high-dimensional text embeddings to interpretable user-defined input label categories. Zero-shot classification reframes the classification task as a language inference problem that determines entailment between a text and the user-supplied label categories yin2019benchmarkingzeroshottextclassification. Notably, this language inference task does not require labeled training data. On the peer effect identification side, we develop new econometric methods based on instrumental variables (IV) to accommodate correlated multivariate outcomes and multidimensional latent variables. Our methods do not require dense networks and accommodate flexible nonparametric modeling of multidimensional latent homophily vectors. The proposed estimator is shown to be $\sqrt{N}$ consistent and asymptotically normal.
We apply this new unified framework of LLM embeddings and econometric methodology to written exchanges between residents of three low-security correctional facilities in a midwestern state in the United States. These correctional facilities, formally known as therapeutic communities (TC), are for individuals with criminal behavior and substance use disorder de2000therapeutic. These programs are based on mutual aid, encouraging residents to be collaborative, practice group therapy, and provide regular feedback to their peers using affirmations and corrections. Corrections are targeted towards individuals for violating TC norms, and affirmations are given to encourage positive and prosocial behavior de2000therapeutic,perfas2014therapeutic. Affirmations and corrections are written on “slips,” which are brief forms containing fields for the sender, receiver, and date, along with the message. Previous randomized controlled trials have demonstrated that TCs reduce criminal activity and drug use sacks2012randomized,bahr2012works.
In these three correctional facilities, we observe the administrative records for each resident's precise entry and exit dates, pre-entry covariates, and details on the content, date, sender, and recipients of these written exchanges. These residents are then mapped with the recidivism records provided by the midwestern state's department of rehabilitation and correction. These correctional facilities are unique both in their practice of exchanging written affirmations and corrections and also in their process of recording them. This unique feature enables us to study peer effects in these language interactions and to use information contained in language to predict meaningful downstream outcomes, such as recidivism.
Using only the pre-entry covariates, we can predict 3-year recidivism with an out-of-sample AUC only slightly above 0.5, which is just marginally better than random guessing. However, using transformer-based LLM embeddings of the 80,000 to 120,000 written exchanges, prediction accuracy improves dramatically to an AUC of 0.70, approximately 30% higher than using pre-entry characteristics alone. This performance is comparable to or slightly better than published AUC values from state-of-the-art commercial software, such as COMPAS, which uses hundreds of individual characteristics in a machine learning model, raising concerns about fairness dressel2018accuracy. An AUC of 0.70 is also generally considered “good” in the context of predicting criminal recidivism laqueur2024algorithmic. These results establish that LLM representations of written exchanges capture economically meaningful behavioral variation. This finding further supports recent work suggesting that peer interactions are key mechanisms behind TC effectiveness campbell2019relationship,nath2022identifying.
After establishing that written exchanges in TCs contain meaningful signals for recidivism, we examine whether these peers influence each other's language profiles. Using zero-shot classification, the unstructured language exchanges are mapped to four user-defined input label categories that represent distinct behavioral dimensions. These chosen input labels are community support, personal growth, rule violations, and disruptive conduct.\footnote{These four input labels for zero-shot classification were chosen using principled, iterative, and replicable LLM prompts; the steps are detailed in Appendix (ref).} Applying the peer effects methodology developed in this paper to the predicted probabilities for language categories reveals significant peer effects. Residents affect each other's language profiles even after accounting for unobserved homophily using the new methods developed in this paper. Peer effects within the same language dimension are substantially larger than cross-category spillovers across dimensions, a pattern observed consistently across both male and female units and across sender and receiver profiles.
The novel econometric methods developed and employed in this study indicate significant peer effects in language profiles, and predictive analysis demonstrates that language exchanges are associated with 3-year recidivism. Nevertheless, predictive power alone is not sufficient for policy implications. Estimating the causal impact of language profiles on recidivism is necessary to determine whether language exchanges in TC can causally affect three-year recidivism. The primary identification challenge in this causal relationship is the presence of unobserved confounders that may influence both language profiles during the stay in TC and recidivism. We estimate the association between language profiles and recidivism, controlling for covariates that we can observe. These covariates include Level of Service Inventory-Revised (LSI-R), which is an aggregate measure of 54 questions provided to participants spanning criminal history, mental health, education, employment, family, living conditions, and alcohol and drug abuse. However, we acknowledge the limitation of this estimate due to unmeasured confounders and use two approaches oster2019unobservable,cinelli2020making to bound the magnitude of omitted-variable bias in the discussion section. We find that estimates of the effect of language in TC on recidivism remain statistically significant at conventional levels, even after accounting for considerable omitted-variable bias relative to the observables.
Text analysis is becoming increasingly important in economics for extracting meaningful insights from unstructured data. The works of kalamara2022making,lippmann2022gender,brehm2025vaccines,alsan2025something use some of the earlier approaches to text analysis, such as sentiment analysis, word frequency counts, bag-of-words, topic modeling, and dictionary methods. Some recent works, such as gennaro2022emotion,ash2022ideas, have improved on these earlier models and instead deployed word embedding methods such as Word2Vec mikolov2013efficient, for text analysis. However, as outlined in the introduction, these pre-transformer word embedding methods do not account for context. In contrast, we deploy transformer-based foundation models to estimate text embeddings and then use them for zero-shot classification, leveraging the attention mechanism of vaswani2017attention. In this regard, our approach is related to bajari2023hedonic, which uses embeddings from transformer models to construct hedonic price indices from Amazon data.
The linear-in-means peer effect model has been widely employed in econometrics for estimating peer effects manski1993identification,bramoulle2009identification. bramoulle2009identification established identification conditions leveraging the network structure and instrumental variables when networks are exogenous conditional on observables. However, networks are typically endogenous due to the phenomenon of latent homophily. These are unobserved characteristics that drive both network formation (peer selection) and the outcomes. Recent work addresses this problem through latent variable models goldsmith2013social,johnsson2021estimation, but existing methods require unrealistically dense networks (expected degree scaling as $N$) and rely on univariate latent homophily. Other work accommodates multivariate correlated outcomes zhu2020multivariate,cohen2018multivariate but assumes exogenous network formation. We extend the peer effects methodology to simultaneously accommodate sparse networks (expected degree scaling as $\sqrt{N}$), multidimensional latent homophily, and multivariate outcomes, proving $\sqrt{N}$-consistency under realistic sparsity conditions and using nonparametric sieve methods to control for latent confounding flexibly.
We also contribute to modeling networks using multidimensional additive and multiplicative latent-variable models for peer-effect estimation. While the works of goldsmith2013social,johnsson2021estimation integrated network formation with peer effects estimation, these approaches have used simplistic network models with univariate additive latent factors to account for the unobserved confounders. We advocate the use of Random dot product graph (RDPG) models athreya2017statistical, which are multiplicative latent variable models, and the latent space models, which accommodate additive and multiplicative latent effects along with observed covariates hoff2002latent,ma2020universal,li2023statistical. In contrast to the one-dimensional latent-variable model of graham2017econometric that johnsson2021estimation employed, the models we deploy use multidimensional latent variables to model network data flexibly. We build methodological and theoretical tools needed to enable the use of this rich class of additive and multiplicative latent space models in peer effect identification, using estimated latent positions with a nonparametric sieve adjustment.
We empirically compare five network specifications (RDPG, latent space with and without covariates, and tetrad logit and joint fixed effect estimators of graham2017econometric) in their ability to replicate observed network properties (modularity, degree heterogeneity, transitivity) of our correctional facility network data. We find that the latent space model with additive and multiplicative factors plus covariates ma2020universal best fits our data.
Our findings have important policy implications for the criminal justice system. Recidivism imposes substantial costs. Reconviction rates range from 18-55% across countries yukhnenko2023criminal, with U.S. incarceration costs exceeding one trillion dollars annually pettus2016economic estimated in 2016. Recidivism prediction has been a central focus in criminal justice policy, with widely-used risk assessment tools such as COMPAS incorporating over 100 pre-entry features dressel2018accuracy. Despite extensive feature engineering, prediction accuracy remains limited dressel2018accuracy,lin2020limits. Our results demonstrate that language exchanges contain rich behavioral signals in correctional settings where pre-entry covariates have extremely limited signals for recidivism.\footnote{While prediction tools inform risk assessment and resource allocation, designing effective interventions requires identifying factors that causally affect recidivism. A separate strand of research has identified several such factors, including employment opportunities, welfare programs, and neighborhood institutions barrios2025recidivism,tuttle2019snapping,galbiati2021jobs.}
This paper uses datasets from two male units and one female unit of Therapeutic Communities (TCs). These TCs in a midwestern state in the United States were minimum-security, community-based correctional facilities. Thus, while they were locked facilities, they were not units in a prison. These TCs had residents from a mixed urban/suburban/rural area; each unit had 80 beds. Since data was collected over time, the female unit had 982 unique residents, and the male units had 1649 unique residents. The two male units are located on two floors of the same TC, and henceforth they will be collectively referred to as the male unit. We observe the entry and exit dates of these residents. The residents may stay in these units for up to 180 days. They may either successfully graduate from the program at the end of their stay or leave it unsuccessfully. We find, on average, the residents stayed for about 120 days. Table (ref) (in the Appendix) details the residents used in the analysis. There is a very high graduation rate from these units, with 87% of residents successfully graduating from the female unit, and the corresponding percentage of graduates for the male unit is 89%.
Given that we observe the time stamps of entry and exit, we observe a substantial variation in these dates and the residents' length of stay in these units. Figure (ref) (a,b) shows the variation in these entry and exit dates, and panel (c) shows substantial differences in the length of the stay. The resident information from TCs is mapped with the administrative dataset maintained by the State's Department of Rehabilitation and Corrections. This administrative dataset provides the date of reincarceration of TC residents, tracked up to 3 years. Reincarceration is one of several possible measures of recidivism, but it is the most commonly used in the criminology literature bouchard2024seeing and has been used in previous studies of TC outcomes doogan2016semantic,warren2007my,warren2020tightly. Table (ref) shows that 22% of the residents had recidivism in the female unit and 25% in the male unit (among all residents, irrespective of graduation/success from TC treatment status). Interestingly, these numbers align with global figures obtained by yukhnenko2023criminal across 33 countries and various types of crimes.
During their stay in TCs, the residents were encouraged to exchange written affirmations and corrections with each other as part of their therapy. We observe the timestamps, text, sender, and receiver IDs for each of these exchanges in these TCs. The digitized records of these texts consist of 123,000 such exchanges between the residents in the female unit and about 83,000 exchanges between the residents in the male unit. Figure (ref) (Appendix) shows the distribution of the count of affirmations (push) and corrections (pull) aggregated by sender ID in panels (a,b) and by receiver ID in panels (c,d), respectively. The female unit shows differences in the distribution of affirmations and corrections, with affirmations being less skewed. However, we see an almost overlapping distribution of corrections and affirmations in the male units.
These corrections and affirmations are short messages exchanged between the residents. Figure (ref) provides a few examples of these messages and the data structure. In this figure, we create a sample of four residents who entered the facility around a similar time and exchanged messages. Note that this figure shows examples of only a few messages exchanged due to space constraints.
We deploy transformer-based LLMs to create text embeddings separately for the 123,000 and 83,000 written exchanges in the female and male units. Precisely, we use the BERT (bert-base-uncased) language model from devlin2019bert. \footnote{The implementation is done using the text package in R kjell2023text that allows users to access these LLMs easily from \href{https://huggingface.co/}{huggingface.co}.} We use these embeddings to create sender and receiver profile vectors for the residents based on the language they use to interact with their peers. These embeddings represent how the behavior and engagement of individuals vary within the TC. The sender profiles indicate how individuals interact with their peers, for example, being supportive or dismissive of others. The receiver profiles, on the other hand, are related to the behavior of the individuals as perceived by their peers. This is because the messages individuals receive reflect their conduct during their stay in the TC (for example, being called out for disruptive conduct). We use these embeddings to assess two key questions- one related to predictive ability for recidivism and the other related to the mechanism of peer influence. We describe in Section (ref) that these word embeddings have substantial predictive ability for recidivism.
However, it is hard to interpret these 768-dimensional embeddings meaningfully. Therefore, our second method uses a LLM-based Zero-Shot approach yin2019benchmarkingzeroshottextclassification,meshkin2024harnessing,wang2023large using the BART model of lewis2019bart. In this method, the user provides a fixed class of words to the LLM along with the text messages, and the method generates predictions for the user-supplied classes for each message. Note that the classification is performed without seeing any labeled training data (and hence “Zero-Shot”), and with only the knowledge of natural language that the LLM has acquired.
The classification mapping is explained in Figure (ref). The figure provides examples of 5 text messages from TCs and their associated LLM-based Zero-Shot classifier probabilities.\footnote{This method is implemented using the transforEmotion package in R christensen2024transforemotion. Similar to the text package, any LLM-based Zero-Shot classifier with a pipeline on Huggingface can be implemented on the local machine without sending any data to an external server.} The four user-supplied classes are selected using LLM prompts to the latest AI models in a principled, iterative, and reproducible manner. See section (ref) in the Appendix for further details on the input label generation. In summary, for choosing the input labels, we generate 100 random samples of text messages stratified by message type. These random samples are then supplied to the Claude API (claude-sonnet-4) with a fixed query for 100 iterations. Finally, the input labels are chosen based on the clusters that well represent the LLM-generated labels in the 100 iterations. Notably, the LLM-based Zero-Shot classification method has been shown to achieve performance comparable to that of deep neural networks with access to labeled training data in biomedical applications meshkin2024harnessing.
We display the scatter plot for class probabilities by message type. The classification appears to work well as the probabilities for “rule violations" and “disruptive conduct" are much higher for corrections, and we see the opposite for affirmations (Figure (ref) for females and Figure (ref) for males, respectively). We review the methodology behind transformers, LLMs, and Zero-Shot classification in Section (ref).
Our final step is to study the underlying behavioral mechanism for these exchanges. Using the low-dimensional interpretable class probabilities, we set up a multivariate outcome model of peer effects that allows for contextual and endogenous peer effects and endogeneity in the network formation. The network links are generated using information on entry and exit date stamps, along with the dates, sender ID, receiver ID, and message exchange frequency. We describe our peer effect estimation methodology in more detail in section (ref).
We describe our methodology that includes five components (described over seven subsections): (1) obtaining embedding vectors of the text exchanges from a pre-trained language model, (2) classifying the messages into user-defined categories through Zero-Shot classification, (3) using the embedding vectors into a high dimensional penalized regression model for predicting recidivism, (4) peer effect estimation using multivariate outcomes and endogenous networks, (5) modeling the network with a sparse mutli-dimensional latent variable model.
We deploy Pre-trained LLMs on the corpus of written messages to derive two representations: (1) context-sensitive embedding vectors and (2) probability vectors for meaningful psychological categories using Zero-Shot classification. We use the bert-base-uncased model devlin2019bert for (1), and the bart-large model lewis2019bart for (2). Before detailing the model structure, we briefly review the development of text embedding methods and discuss our rationale for adopting transformer-based AI models, such as BERT and BART.
A text embedding is a high-dimensional numerical vector representation of text. If two words convey similar meanings, they are represented closely in the vector space. However, earlier word embedding models, such as Word2Vec mikolov2013efficient, did not consider the context for words, limiting their ability to handle words with multiple context-specific meanings. Transformer-based models, like BERT and BART, overcome these limitations using self-attention mechanisms vaswani2017attention, where the embedding for each word is dynamically influenced by the embedding of the entire input sequence, providing potentially different embeddings for each word based on context. The generative pre-training (GPT) framework Radford2018ImprovingLU uses semi-supervised learning in the transformer architecture. The model is pre-trained unsupervised on a large data corpus to learn embeddings or representations, followed by supervised fine-tuning for specific tasks. The language representation model BERT in devlin2019bert further improved the capabilities of pre-trained language models by introducing a bi-directional transformer, which can use both left and right context for a word as opposed to unidirectional processing (typically from left to right) in OpenAI GPT-1 of Radford2018ImprovingLU.
The first step to processing text is tokenizing it into words or sub-words $z_k$ for $k = 1, \ldots, K$, forming a sequence $\mathbf{z} = (z_1, \ldots, z_K)$ of tokens where $K$ is the total number of tokens for each written text. For each token position $k$, both BART and BERT construct an input embedding as a sum of token embedding $\mathbf{e}_k^{\mathrm{token}}$ and a positional embedding $\mathbf{e}_k^{\mathrm{pos}}$ as $\mathbf{e}_k = \mathbf{e}_k^{\mathrm{token}} + \mathbf{e}_k^{\mathrm{pos}}.$ In BERT, an additional segment embedding $\mathbf{e}_k^{\mathrm{seg}}$ is included if input text consists of a sentence pair: $\mathbf{e}_k = \mathbf{e}_k^{\mathrm{token}} + \mathbf{e}_k^{\mathrm{pos}} + \mathbf{e}_k^{\mathrm{seg}}.$ These vectors then form an input embedding matrix $\mathbf{E} = [\mathbf{e}_1, \ldots, \mathbf{e}_K]^\top \in \mathbb{R}^{K \times H}$. Here $H$ denotes the dimension of each embedding vector. Precisely, $H = 768$ for bert-base-uncased and $H = 1024$ for bart-large.
A sequence of transformer blocks processes $\mathbf{E}$, where the output of each block is fed as an input into the next block. Let the output of the $l$-th block be $\mathbf{E}^{(l)} \in \mathbb{R}^{K \times H}$, with the initial input being $\mathbf{E}^{(0)} = \mathbf{E}$. We review the processes involved in each of these transformer blocks and the multi-head self-attention mechanism in section (ref) of the Appendix.
To numerically encode the semantic content of each message, we use the output from the 11th encoder block (i.e., the penultimate block) of bert-base-uncased model to construct an “embedding profile" for each participant --- a high-dimensional vector summarizing their language use across multiple messages during their stay in TC. We choose the penultimate layer as that layer is widely thought to contain the unsupervised representation of data in the representation learning framework bengio2013representation, while the last layer is trained or fine-tuned to be task-specific. One can, in principle, also choose an embedding representation from an earlier layer. Let $\mathbf{f}^{\mathrm{token}}_{c,k} \in \mathbb{R}^H$ denote the embedding of the $k$-th token in the $c$-th message, extracted from the 11th encoder. We obtain the message-level embedding by taking the average of these vectors over all tokens in the message as $\mathbf{f}^{\mathrm{msg}}_c=\frac{1}{N_c} \sum_{k=1}^{N_c} \mathbf{f}^{\mathrm{token}}_{c,k}$, where $N_c$ is the number of tokens in message $c$. We then compute the “sender embedding profile” for a sender $s$ by taking an average of the the message-level embeddings, $\mathbf{f}^{\mathrm{sender}}_s=\frac{1}{N_s} \sum_{c \in \mathcal{C}_s} \mathbf{f}_c^{\mathrm{msg}}$, where $\mathcal{C}_s$ denotes the set of messages sent by sender $s$, and $N_s = |\mathcal{C}_s|$ is the number of such messages. The embedding profile for each receiver can be obtained analogously.
To interpret the embeddings, we relate them to interpretable psychological constructs, using the “Zero-Shot" classification framework yin2019benchmarkingzeroshottextclassification. In Zero-Shot classification, a user provides an arbitrary number of labels, and the LLM assigns given texts to the most appropriate label. This enables us to assign written affirmations and corrections to psychologically relevant categories without requiring additional training.
In Zero-shot learning, the classification problem is turned into a language inference problem. We use bart-large model fine-tuned on the MNLI (Multi-Genre Natural Language Inference) task dataset, a dataset designed for language inference task. In MNLI, each input consists of a premise and a hypothesis, and the pre-trained BART model is trained to classify the semantic relationship between them as one of three predefined MNLI classes: entailment, neutral, or contradiction. A premise-hypothesis combination is concatenated into a single input sequence using special tokens following the format: <s> premise </s> hypothesis </s>, which allows the model to process both segments jointly while preserving their roles.
The encoder uses a self-attention mechanism that allows the model to capture contextual dependencies between tokens across both segments. This enables the model to assess how strongly the premise supports, contradicts, or is unrelated to the hypothesis. During training, the embedding of the first token <s> from the final encoder layer is passed to a classification head, which outputs a logit vector over the three MNLI classes.
Figure (ref) illustrates how a TC text message is classified into user-defined categories using this framework. In our application, each peer message is treated as a premise, and we define four psychologically meaningful word strings: “Rule Violations", “Disruptive Conduct", “Community Support", and “Personal Growth". Each of the four strings is converted into a hypothesis using a fixed template: “This example is \{\}."
For each message $c$, the model generates four premise-hypothesis pairs by combining the message with each label-based-hypothesis indexed by $\ell=1, \ldots,4$. For each pair, the model obtains a probability vector for entailment $\mathbf{s}_{c} =(p_{c,\ell}^{\mathrm{entail}})_{\ell=1}^4 \in \mathbb{R}^4$, which reflect how strongly the input message supports each of the user-defined categories. These four entailment scores are then passed through a softmax function to obtain a probability vector over the labels: $\boldsymbol{\pi}_c = \mathrm{softmax}(\mathbf{s}_c)$. Finally, we aggregate these scores at the level of individuals to create personal class profiles, analogous to the embedding-based profiles described earlier, but in a more tractable lower-dimensional space.
We fit a penalized regression lasso method tibshirani1996regression to predict recidivism using the resident-level covariates along with the individual embedding profiles. Suppose $Y_i$ denotes the recidivism status for individual $i$, which takes the value $1$ if the individual recidivates and $0$ if the person does not recidivate. Therefore, we can model $Y_i$ with a Bernoulli distribution with probability of success being $p_i$ and then further model $p_i$ using the covariates and embedding profiles. However, since the embedding profiles are high-dimensional, we use a penalized logistic regression method (logistic-lasso). Suppose $\bm X$ is the $N\times p$ matrix of observed covariates and $\bm T$ is the $N\times H$ matrix containing the dimensions of text embeddings obtained from the LLM model as columns. Then the logistic regression model is $\log \left(\frac{p_i}{1-p_i}\right)=\beta_{0}+\bm X_i\beta_{1}+\bm T_i\beta_{2} $, where $\beta_1$ and $\beta_2$ are $p$ and $H$ dimensional parameter vectors respectively. We optimize for $\beta_1, \beta_2$ by maximizing the penalized likelihood function with $\ell_1$ penalty as follows, \[ (\hat{\beta}_1, \hat{\beta}_2) = \text{argmin} \{ - \ell (\beta_1,\beta_2, \bm X, \bm T) + \lambda_1 \|\beta_1\|_1 + \lambda_2\|\beta_2\|_1\}. \]
The parameters $\lambda_1, \lambda_2$ denote the penalty parameters which maybe different for the observed covariates and the embedding vectors. We assess the model's predictive accuracy with out-of-sample AUC values through a five-fold cross-fitting.
Beyond predicting recidivism, we also want to interpret the link between individual behavioral and emotional profiles as manifested by the messages with recidivism. In order to do so we perform a Zero-Shot classification as described before to classify the messages into the four groups. Let $\bm Q$ denote the $N \times 4$ matrix containing the average of probabilities assigned to the messages sent (respectively received) by each individual. Since the predictors are now low-dimensional (just four dimensions in addition to the covariates), we will not need to penalize the coefficients. However, we note that some of the covariates are now “compositional", since $\sum_j Q_{ij}=1$ for each resident $i$. Therefore, following standard techniques of dealing with compositional predictors, we designate one covariate as baseline and take log ratios (ALR) of other covariates with respect to the baseline covariate aitchison1982statistical. This is important for interpreting the model coefficients, but not for prediction accuracy. \footnote{We note that we could just continue with the original predictors and drop the intercept, and have a model with comparable prediction accuracy.}
We consider the multivariate peer effect model, which captures dependencies among multiple outcomes observed across spatially or network-connected entities. Suppose we have $N$ entities and we observe $m$ outcome variables, and $p$ covariates for each entity. Further, we observe a network among the entities whose adjacency matrix is $\mathbf{A}$, where the elements $a_{ij}$ represent the relationship between entities $i$ and $j$. Then we define the multivariate peer effect model as:
where $\mathbf{Y} \in \mathbb{R}^{N \times m}$ is the matrix of outcome variables, $\mathbf{G}=(g_{ij})\in\mathbb{R}^{N \times N}$ where $g_{ij}:=\frac{a_{ij}}{\sum_{j\ne i}a_{ij}}$, is the row-normalized counterpart of the network adjacency matrix $\mathbf{A}$, and $\mathbf{X} \in \mathbb{R}^{N \times p}$ is the matrix of observed covariates. We write $\bm Y_i\in\mathbb{R}^m$ and $\bm Y_{.,j}\in\mathbb{R}^N$ to denote the $i$-th row and $j$-column vectors of an $N \times m$ matrix $\mathbf{Y}$, respectively. Matrix $\mathbf{D} \in \mathbb{R}^{m \times m}$ is the peer effect parameter matrix. The diagonal of $\mathbf{D}$ are the direct peer spillover effects on the same dimensions of $\mathbf{Y}$, and the off-diagonal elements are the indirect peer spillover effects of a dimension on another dimension. The matrices $\mathbf{B}_1, \mathbf{B}_2 \in \mathbb{R}^{p \times m}$ are coefficient matrices for $\mathbf{X}$ and $\mathbf{G}\mathbf{X}$, respectively.
However, since the network is endogenous, we cannot identify and consistently estimate $\mathbf{D}$ from this model using bramoulle2009identification's IV 2LSLS framework. For concreteness, we assume for each individual there is a vector $\mathbf{U}_i$ of latent homophily variables that is correlated with both the network adjacency matrix $\mathbf{A}$ and the error matrix $\mathbf{E}$. We assume the error matrix $\mathbf{E} \in \mathbb{R}^{N \times m}$ is such that for each row $\mathbb{E} [ \mathbf{E}_i| \mathbf{U}_i] =\bm{h}^E(\bm U_i)$ for some unknown function $\bm h^E$, and $ \mathbb{V}[ \mathbf{E}_i| \mathbf{U}_i] = \mathbf{V}$. The matrix $\mathbf{V}$ allows dependence across outcome dimensions. The vectors $\mathbf{E}_i,\mathbf{X}_i,\mathbf{U}_i$ are assumed to be i.i.d. across individuals $i$. We let $A_{ij}=f(\mathbf{U}_i,\mathbf{U}_j,\xi_{ij})$, where $\xi_{ij}$s over all $(i,j)$ are assumed to be i.i.d and independent of $\mathbf{X,E,U}$. We postpone the discussion on the function $f$ and the models for network formation until the next section. We further assume that $\mathbb{E} [\bm E_i | \bm X_i, \bm U_i] = \mathbb{E}[\bm E_i | \bm U_i]$, i.e., conditional on $\bm U_i$, the covariates $\bm X_i$ and error term $\bm E_i$ are uncorrelated.
This model is similar to the multivariate extension of the spatial autoregressive (MSAR) model proposed in zhu2020multivariate. However, the model in zhu2020multivariate can be used for peer effect estimation only under the assumption of exogenous network formation ($\mathbf{G}$ is uncorrelated with $\mathbf{E}$) and the proposed estimator is the maximum likelihood estimator and not IV2SLS. Further, model ((ref)) can also be thought of as extending the endogenous and simultaneous peer effect model in johnsson2021estimation, both in terms of multidimensional outcome and latent variables.
To identify the parameter matrices, we take an instrumental variable approach similar to johnsson2021estimation, but extend the methodology in several directions as we detail below. We define the combined regressor and instrument matrices as: \[ \mathbf{Z} = [\mathbf{G}\mathbf{Y}, \mathbf{X}, \mathbf{G}\mathbf{X}] \in \mathbb{R}^{N \times (m + 2p)}, \quad \mathbf{K} = [\mathbf{X}, \mathbf{G}\mathbf{X}, \mathbf{G}^2\mathbf{X}] \in \mathbb{R}^{N \times 3p}. \] The stacked parameter matrix is: \[ \boldsymbol{\beta} = [\mathbf{D}, \mathbf{B}_1, \mathbf{B}_2]^T \in \mathbb{R}^{(m + 2p) \times m}. \] Using this notation, equation (ref) can be written compactly as $ \mathbf{Y} = \mathbf{Z} \boldsymbol{\beta} + \mathbf{E}.$ A natural estimator for $\boldsymbol{\beta}$ is the two-stage least squares (2SLS) estimator:
This estimator is a multivariate extension of the estimator proposed in bramoulle2009identification,kelejian1998generalized and is a consistent estimator provided the network is exogenous, which is not the case in our setup. Taking expectations on both sides of Equation (ref) conditioning on $\bm U$ and subtracting this from equation (ref), we obtain
If we redefine, $\tilde{\bm Y} = \bm Y- \mathbb{E}[\bm Y|\bm U]$, $\tilde{\bm Z} = \bm Z-\mathbb{E}[\bm Z|\bm U]$, and $\tilde{\bm K} = \bm K - \mathbb{E}[\bm K| \bm U].$ Then we can write the above equation as, $ \tilde{\mathbf{Y}} = \tilde{\mathbf{Z}} \boldsymbol{\beta} + \tilde{\mathbf{E}},$ where $\tilde{\mathbf{E}} = \mathbf{E} - \mathbb{E}[\mathbf{E} \mid \mathbf{U}]$. For this redefined linear model we can solve for $\bm \beta$ using IV-2SLS with the redefined instrument matrix $\tilde{\mathbf{K}}$. The moment equation for this is given by
The following proposition provides sufficient conditions for the identification of the true parameter matrix $\beta_0$.
Therefore, given access to the latent positions $\mathbf{U}$ and the three conditional mean functions, the parameter $\boldsymbol{\beta}_0$ can be estimated using 2SLS. Note that the rank condition in Proposition (ref) can generally be satisfied with $p>m$, i.e., if we have more covariates than dimensions. However, the moment equation still cannot be estimated as we do not observe $\mathbf{U}$ and do not know the conditional expectations of $\mathbf{Y,Z,K}$ given $\mathbf{U}$. We next address how to estimate both of these components.
As mentioned earlier, we model the network data using a latent variable model $A_{ij}=f(\mathbf{U}_i,\mathbf{U}_j,\xi_{ij})$ and posit that the node level latent variables $\mathbf{U}_i$ are involved in both the network formation model and are correlated with the error term of the outcome model. These latent variables represent unobserved characteristics of individuals, which may be responsible for tie formation (latent homophily), and an unknown function of these latent variables is part of the error term in the outcome model. Therefore, these latent variables can be estimated from the observed network. This is a core assumption made in recent methodological advances in peer effect estimation mcfowland2021estimating,johnsson2021estimation,nath2022identifying,goldsmith2013social. However, while mcfowland2021estimating,nath2022identifying focused on the longitudinal peer effect model, we consider the simultaneous peer effect model. Further, mcfowland2021estimating,nath2022identifying assumed the latent variables enter the model for $\mathbf{Y}$ linearly as $\mathbf{U\beta}$. In contrast, we assume the conditional expectations of $\mathbf{Y,Z,K}$ to be unknown functions $\mathbf{h(U)}$ of the latent variables. In this aspect, our approach resembles johnsson2021estimation; however, we consider more general latent variable models for network data with multi-dimensional latent variables, which are more realistic for modeling individual characteristics. Further, the theoretical results in johnsson2021estimation require the network to be “dense” in the sense that every individual is expected to be connected to $O(N)$ individuals and consequently, the total number of edges in the network is expected to be $O(N^2)$. This is quite unrealistic as real-world networks tend to be sparse with individuals connected to only a few other people, even when the network size is huge (e.g., in online social networks, even if there are millions of people present in the network, an individual member only has a few hundred or thousands of connections). In contrast, our results here hold for sparse networks.
We define the stacked data matrix $\mathbf{W} := [\mathbf{Y}, \mathbf{Z}, \mathbf{K}] \in \mathbb{R}^{N \times (2m+5p)}$. The conditional mean functions are defined as $\bm{h}^Y(\bm{U}) := \mathbb{E}[\bm{Y} \mid \bm{U}]$, $\bm{h}^Z(\bm{U}) := \mathbb{E}[\bm{Z} \mid \bm{U}]$, and $\bm{h}^K(\bm{U}) := \mathbb{E}[\bm{K} \mid \bm{U}]$. These are collected into the full conditional mean matrix $\bm{h}(\bm{U}) := [\bm{h}^Y(\bm{U}), \bm{h}^Z(\bm{U}), \bm{h}^K(\bm{U})] \in \mathbb{R}^{N \times (2m+5p)}$. The function $\mathbf{h}: \mathbb{R}^d \to \mathbb{R}^{2m+3p}$ is an unknown non-linear function of $\mathbf{U}_i$. Let \( h_j(\cdot):\mathbb{R}^N\to\mathbb{R} \) denote the \( j \)th component of \( \bm{h}(\cdot) \), for \( j = 1, \dots, (2m+5p) \). Then, for each \( \bm{U}_i \), \( h_j(\bm{U}_i) \in \mathbb{R} \) is the scalar value of the \( j \)th conditional mean component. In what follows, to keep the notation simple, let $\mathbf{u}$ denote a generic row $\mathbf{U}_i$.
To estimate each component function $h_j(\bm{u})$, we use tensor-product basis functions zhang2023regression given by $\phi_k(\bm{u}) = \prod_{\ell=1}^{d} \psi_{k_\ell}^{(\ell)}(u_\ell)$, where each $\psi_{k_\ell}^{(\ell)}$ is a univariate basis function (polynomial or cosine) on the $\ell$-th coordinate of $\mathbf{u}$. With these basis functions, we approximate each component function \( h_j(\bm{u}) \), by a linear combination of basis functions: \[ h_j(\bm{u}) \approx \sum_{k=1}^{L_N} \phi_k(\bm{u}) \, \alpha_k^j = \sum_{k=1}^{L_N} \prod_{\ell=1}^{d} \psi_{k_\ell}^{(\ell)}(u_\ell) \, \alpha_k^j , \] where \( \bm{\alpha}^j = (\alpha_1^j, \dots, \alpha_{L_N}^j)^\top \in \mathbb{R}^{L_N} \) is the coefficient vector, and the number of basis functions to use $L_N$ is possibly a function of $N$. Note this formulation allows for different univariate basis functions indexed by $\psi_{k_\ell}$ for each component basis $\phi_k$ and each dimension of $u_l$. It also allows for interactions among the predictors.
For the generic vector $\mathbf{u}$, let \( \bm{\phi}^{L_N}(\bm{u}) := \left( \phi_1(\bm{u}), \dots, \phi_{L_N}(\bm{u}) \right)^\top \in \mathbb{R}^{L_N} \). Then we construct the design matrix for all $N$ observations as: $\mathbf{\Phi}_N := \left[\bm{\phi}^{L_N}(\bm{u}_1), \ldots, \bm{\phi}^{L_N}(\bm{u}_N)\right]^T \in \mathbb{R}^{N \times L_N}$. The coefficient vector \( \bm{\alpha}^j \) is obtained by regressing the observed vector \( \mathbf{w}_{\cdot j} \) on the basis matrix \( \mathbf{\Phi}_N \) via ordinary least squares, as $\hat{\bm{\alpha}}^j = (\mathbf{\Phi}_N^\top \mathbf{\Phi}_N)^{-1} \mathbf{\Phi}_N^\top \mathbf{w}_{\cdot j}$. The fitted values for the \( j \)th component function at all sample points are then given by: $ \hat{h}_j(\bm{U}) := \mathbf{\Phi}_N \hat{\bm{\alpha}}^j = \mathbf{P}_{\mathbf{\Phi}_N} \, \mathbf{w}_{\cdot j},$ where \( \mathbf{P}_{\mathbf{\Phi}_N} := \mathbf{\Phi}_N (\mathbf{\Phi}_N^\top \mathbf{\Phi}_N)^{-} \mathbf{\Phi}_N^\top \) is the projection matrix associated with the basis space, and \( A^{-} \) denotes any symmetric generalized inverse of $A$.
Since latent positions $\bm{U}$ are unknown, we use a two-step procedure. For estimating $\bm U $, we model the network using latent variable models such as the RDPG model athreya2017statistical,rubin2022statistical,xie2023efficient and the additive and multiplicative effects latent space model ma2020universal,hoff2021additive,hoff2002latent,li2023statistical. Then our two-stage procedure is as follows. (i) we first construct an estimator \( \widehat{\bm{U}} := (\widehat{\bm{u}}_1, \dots, \widehat{\bm{u}}_N)^\top \) for the latent traits using spectral embedding for RDPG models or maximum likelihood estimation for additive and multiplicative effects latent space models, and (ii) we evaluate the basis functions at \( \widehat{\bm{U}} \) and apply the same projection strategy. With these estimated latent vectors, we define the estimated design matrix as $ \widehat{\mathbf{\Phi}}_N := \mathbf{\Phi}_N(\widehat{\bm{U}}) = \left[ \bm{\phi}^{L_N}(\widehat{\bm{u}}_1), \ldots, \bm{\phi}^{L_N}(\widehat{\bm{u}}_N) \right] \in \mathbb{R}^{N \times L_N}.$ The final estimated conditional mean functions are $\hat{\bm h}^W(\widehat{\bm{U}}) := \mathbf{P}_{\widehat{\mathbf{\Phi}}_N}\bm W$ where $\mathbf{P}_{\widehat{\mathbf{\Phi}}_N} := \widehat{\mathbf{\Phi}}_N (\widehat{\mathbf{\Phi}}_N^\top \widehat{\mathbf{\Phi}}_N)^{-1} \widehat{\mathbf{\Phi}}_N^\top$.
Denote $\bm M_{\hat{\bm \Phi}_N} = \bm I_N - \bm P_{\hat{\mathbf{\Phi}}_N}$. Then, our Two-Stage Least Squares (2SLS) estimator is:
In terms of latent variable network models, we first consider the Latent Space Model (LSM), which is a general random graph model that includes both additive and multiplicative latent variables and can also accommodate covariates hoff2002latent,hoff2021additive,ma2020universal,li2023statistical. The model is parameterized by $d$-dimensional unknown vectors $\bm q_i$ and unknown scalars $v_i$ for all nodes $i=1, \ldots, n$.
We also assume that we have edge level covariates $\bm X= (x_{ij})$ available to us. A network adjacency matrix $\bm A= (a_{ij})$ from this model is generated as follows li2023statistical,ma2020universal: \[a_{ij}\overset{ind.}{\sim}f_{ij}:= f(a;\theta_{ij}),\qquad \theta_{ij}:=\sigma({\bm q_i}'\bm q_j+v_i+v_j + x_{ij}\beta), \quad 1\leq i<j\leq n\] where $f$ is a family of distributions satisfying certain smoothness conditions, and $\sigma:\mathbb{R}\to\mathbb{R}$ is a known link function. In this paper, we consider a particular sparse version of the model due to li2023statistical that assumes $f$ to be a Bernoulli distribution, link function $\sigma$ to be the logistic function, and there exists a sparsity parameter $\rho_N$ such that $\theta_{ij}=logistic({\bm q_i}'\bm q_j+v_i+v_j + x_{ij}\beta +\rho_N)$. Defining $\omega_N = \exp(\rho_N)$, we assume $\omega_N \to 0$ and $\omega_N=\omega(N^{-1/2})$. We let $\bm u_i\triangleq(\bm q_i',v_i)'$ denote the $d+1$ dimensional parameter vector containing all latent variables associated with node $i$. A special case of this model is Random Dot Product Graph (RDPG) model, which can be expressed as follows, \[(a_{ij}|\bm u_i,\bm u_j)\overset{ind.}{\sim}\operatorname{Bernoulli}(\rho_N \bm u_i'\bm u_j),\] where $\rho_N$ again controls the sparsity of the model. We will assume $\rho_N =\omega(\frac{\log^4N}{N})$. Note the expected density (i.e., probability of an edge) scales with $N$ as $\rho_N$, or in other words, the expected degree scales as $N\rho_N$ in this model, making the model suitable for sparse networks. This model does not account for the additive latent variables $v_i$'s but contains $d$-dimensional multiplicative latent variables, and the link function $\sigma(\cdot)$ is specified as identity. Clearly, both these models allow for multidimensional latent variables in the network formation model and are more general than the single latent variable model considered in johnsson2021estimation. Both models are also capable of modeling sparse networks, which the models in johnsson2021estimation,graham2017econometric,auerbach2022identification are not capable of. We compare the model fit of the latent-space network models with the network model considered in johnsson2021estimation,graham2017econometric for the network data in our real-data analysis application (see section (ref)).
We employ the maximum likelihood estimator with Lagrange adjustment in li2023statistical for the estimation of the LSM model. Let $\hat{\bm U}$ denote the maximum likelihood estimator over the constrained parameter space \[\Xi_d:=\{\bm U\in\mathbb{R}^{N\times(d+1)} \,;\,{\bm Q}'\bm 1_N = \bm 0_d, \,{\bm Q}'{\bm Q} \text{ is diagonal, } \|\bm U\|_{2,\infty} = O(1)\text{ as }N\rightarrow \infty\}.\]
Under the set of assumptions laid out in li2023statistical, an upper bound on the error rate of estimating the latent positions in Frobenius norm li2023statistical is $\frac{1}{N}\|\hat{\bm U}-\bm U\|_F^2=O_p(\frac{1}{N\omega_N})$.
Although the LSM described above contains the RDPG model as its special case, there exists a rich amount of results with better convergence rates of the estimate of $\bm U$ when one uses spectral embedding to estimate $\bm U$ in the case of RDPG model athreya2017statistical,cape2019signal,xie2023efficient,chang2024embedding. Specifically, the spectral embedding of the adjacency matrix $A$ can be written as $\hat{\bm U} = \bm U_A|\bm S_A|^{1/2}$, where $\bm{U}_A$ denotes the matrix containing the leading $d$ eigenvectors of $\bm{A}$, which are assumed, without loss of generality, to be ordered by decreasing absolute eigenvalue magnitude. The matrix $|\bm{S}_A|$ is a diagonal matrix containing the absolute values of the corresponding eigenvalues. Then, a faster convergence rate than LSM can be obtained in terms of the $2\to \infty$ norm (maximum of row-wise $\ell_2$ norms) as $\|\hat {\bm U} - \bm U\bm H\|_{2,\infty} = O_{hp}\left(\frac{\log^c N}{N^{1/2}}\right)$ rubin2022statistical,cape2019signal. Here we mean $X_n=O_{hp}(1)$ by: for every $c>0$ there exist constants $M(c),n_0(c)>0$ such that $\Pr(|X_n|>M)<n^{-c}$ for all $n>n_0$.
Our method, combining all the above components, is described in Algorithm (ref). Note that the step estimating $\hat{U}$ with spectral embedding in the algorithm will be replaced with MLE when the latent space model is used for the underlying network.
For a matrix $A$ we denote its spectral norm as $\|A\|=\max_{\|x\|=1}|Ax|$, Frobenius norm as $\|A\|_F=\sqrt{\sum_{ij}a_{ij}^2}$ and the $ 2 \to \infty$ norm as $\|A\|_{2, \infty}=\max_{i}\sqrt{\sum_j a_{ij}^2}$. We use the stochastic order notation $o_p(1)$ to mean that if $x_N = o_p(1)$, then $x_N \overset{p}{\to} 0$. We start with a list of assumptions. The first assumption is a collection of conditions described above in the description of the model, making the model well-behaved and identification possible, similar to those in johnsson2021estimation.
We note again that the network $\bm A$ is endogenous by the dependence between $\bm U_i$ and $\bm E_i$. A consequence of Assumption (ref) is that
i.e., covariates and observed network are uncorrelated with the model error conditional on the latent variables $\bm U_i$.
The above assumption is an extension of the assumption in johnsson2021estimation to multivariate basis functions which says the true component-wise functions $h_j(u)$ are well approximated by the sieve approximations for all vectors $u$ and all components $j$.
Our last assumption restricts $\{\phi_k;k\le L_N\}$ to be a $\zeta_1(k)$-Lipschitz where the order of $\zeta_1(1),\dots,\zeta_1(L_N)$ is determined by the random graph model and appropriate sparsity parameters we assume.
We derive the asymptotic distribution of the 2SLS estimator \( \hat{\boldsymbol{\beta}}_{v,\text{2SLS}} \). Here, we write $\delta_{\max}:=\max_{i\le N}\sum_{j\ne i}a_{ij}$ and $\delta_{\min}:=\min_{i\le N}\sum_{j\ne i}a_{ij}$ to denote the maximum and minimum node degrees, respectively.
The proofs of all results are in the Appendix (ref). Several additional technical lemmas are needed, which are all stated and proved in the Appendix (ref).
Note that the conditions on the graph density indicate that consistent and asymptotically normal estimation is possible even for sparse graphs. Under a binary graph ($a_{ij}$s are either 1 or 0), for example, the first condition in (ref) only requires $\delta_{\max}/\delta_{\min}^2=o_p(1/\sqrt{N})$. A sufficient condition for the second part in (ref) can be (since for a binary graph $\|\bm a_{i_1}\|^2=\delta_{i_1}$ for any $i_1$) \[\mathbb{E}\left[\frac{(\delta_{i_3}\delta_{i_4})^{1/2}\delta_{\max}}{\delta_{i_1}\delta_{i_2}\delta_{\min}^4}\right]\le \mathbb{E} \left(\frac{\delta_{\max}^2}{\delta_{\min}^6}\right)=o(1/N^2).\] These conditions hold under a moderate expected density of the graph. For example, if we let $\delta_{\min} \asymp \delta_{\max}$ and the average expected density of the graph grows as $O(\frac{1}{\sqrt{N}})$, then all conditions are satisfied. This is in contrast to the results in johnsson2021estimation, which only hold if the average expected density of the graph grows as $O(1)$, making the graph unrealistically dense.
An intermediate lemma en route to proving the main theorem further clarifies the distinction between the density requirements for the RDPG and the LSM models.
Clearly, Lemma (ref) provides a result on the concentration of the sieve design matrix when the true latent positions $\bm U$ are replaced with their estimated counterparts $\hat{\bm U}$. Both the RDPG and the LSM models can accommodate sparse networks through the sparsity parameters $\rho_N$ and $\omega_N$ respectively. However, our result with the RDPG model has a better concentration of the $\hat{\bm \Phi}_N$ to $\bm \Phi_N$, for whenever $\omega_N$ is smaller than $\frac{1}{\log^{2c}N}$, i.e., if the network is even slightly sparse. However, we emphasize that our results for both models accommodate sparse networks with density requirement for RDPG being $\frac{\log^4N}{N}$ and LSM being $\frac{N^{1/2}}{N}$ (in addition to what is required to satisfy conditions (ref) amd (ref) in the main theorem).
In this section, we describe a simulation study we performed to evaluate the finite-sample performance of the proposed estimation strategy in the presence of latent homophily and network endogeneity. Specifically, we consider two main scenarios for network and outcome generation: (i) a Random Dot Product Graph (RDPG) model and (ii) a latent space model with covariates. We generate network data, observed covariates, latent variables, and outcomes for each scenario according to the corresponding data-generating process. The simulated data are then used to compare several alternative adjustment strategies for latent confounding within the 2SLS estimation framework. We evaluate both nonparametric sieve-adjustment methods, which use polynomial or tensor product polynomial basis expansions, as well as linear functions of the estimated latent positions. As our primary evaluation metric, we use the mean squared error, calculated across simulation replications.
Table (ref) documents the performance of the estimator when the network is generated from the RDPG model and the spectral embedding is used to estimate the latent variables. The dimension of latent variables in the outcome model is fixed at 2. There are two correlated outcomes where the data generation follows Equation (ref) with endogeneity in network formation. The peer influence parameter matrix is set at $
$ where the diagonal elements capture the direct peer influence spillovers and the off-diagonal elements capture the indirect spillovers in the outcome. The dimension of the coefficients on $h^Y(U)$ varies across the three scenarios. For each scenario, the data is generated 100 times, followed by estimation of $\beta_{v}$ using the steps in algorithm (ref). The table shows that the mean squared error steadily declines with the number of nodes in the network. This supports the fact that the estimator works well in recovering the true parameter of interest as the sample size increases.
Next, we change the data generation and estimation of latent variables, where latent space models that allow for both multiplicative and additive unobserved variables, along with observed covariates, are used for generating the network. A Projected Gradient Descent (PGD) algorithm with Universal Singular Value Thresholding (USVT) initialization is used to estimate the MLE $\hat{U}$ ma2020universal. Note that even though this solves a non-convex optimization problem, the solution has nice theoretical guarantees which we have used in the proofs above. The covariate that explains the network formation is not part of the data generation for the outcomes. The network formation consists of one additive latent factor, two multiplicative latent factors, and two observed covariates. The outcome data generation consists of $[AY,X,AX,h^{Y}(U)]$ where $AY$ is $N\times 2$ matrix, $X$ is $N\times 5$, $AX$ is $N\times 5$ and dimension of $h^{Y}(U)$ varies across the 3 scenarios. The peer influence matrix $D$ used in the data generation for the latent space model with covariates is fixed at $
$. Table \ref{simsfinallsmcov} reports the mean squared error (MSE) across 50 simulations for each scenario as the sample size increases from 100 to 300 to 500. It shows that the MSE consistently declines as the sample size increases for all four parameters of the peer influence $D$ matrix for all three scenarios.
We first investigate the ability of text exchanges to predict recidivism. Recall, we estimate the LLM embeddings for the messages and use those embeddings as predictors in a high-dimensional regression LASSO prediction method. LASSO is implemented with cross-fitting to do the entire prediction exercise out-of-sample.\footnote{Cross-fitting divides the sample into K-folds, and the prediction for the residents in the $k^{th}$ fold uses the model trained on the data from the remaining $ K-1 $ folds.} The penalty parameter is chosen using cross-validation on the entire dataset, and then it is saved and passed into the estimation and prediction for each fold in cross-fitting.
Figure (ref) provides the Receiver Operator Characteristic (ROC) curves and the Area Under the Curve (AUC) for two models, one with both pre-entry covariates ($X$) and embeddings ($T$), and an alternative model without the text embeddings (only $X$) as predictors. Panel (a) shows that the prediction AUC is 30% higher for recidivism in the female unit when we incorporate the LLM-based text embeddings relative to the model with only pre-entry covariates $X$. The dimension reduction is substantial for the former model, as about 26-37 predictors out of text embeddings (768) and pre-entry covariates remain across the five folds in LASSO model fit. Panel (b) compares the prediction accuracy for the male units with and without the embeddings. The prediction accuracy improvement for the males is 30.6%. Figure (ref) repeats the same analysis but for the case when the text embeddings are aggregated by the receiver ID. We find comparable improvements (32.6% for females and 32% for males) in predictive accuracy when the text embeddings are aggregated by the receiver profiles and included as predictors in the LASSO model along with the pre-entry covariates. Recall from section (ref) that the sender profile indicates how an individual interacts with their peers, while the receiver profiles indicate how others perceive an individual.
The out-of-sample AUC values from our model with embeddings are in the range of 0.67-0.70, which is generally considered a good predictive accuracy in the context of recidivism prediction laqueur2024algorithmic. We further note that our pre-entry covariates include an aggregate measurement known as Level of Service Inventory-Revised (LSI-R), a widely used measure in the context of criminal psychology andrews2014psychology. The LSI-R survey contains 54 questions pertaining to criminal history, education, employment, family, alcohol and drug problems, living conditions, and emotional health. Several questions that are part of the LSI-R are known to be predictors of recidivism dressel2018accuracy. Yet, in our sample, the pre-entry covariates, including LSI-R, were only mildly predictive of recidivism. This could be because LSI-R is an aggregated measure, and we do not observe the component scores in our data. Nevertheless, this underscores the importance of AI-based text measurements as potential predictors of recidivism, especially in low-security correctional facilities, where observed covariates have limited ability to predict recidivism.
Next, prediction accuracy is assessed using the class probabilities generated by the transformer-based Zero-Shot classifier. This analysis uses the classes “Personal Growth”, “Community Support”, “Rule Violations”, and “Disruptive Conduct”. Each of these classes is chosen to be indicative of behaviors that are either meaningfully detrimental to community peace and cohesion or supportive of personal growth and community well-being. We use the class probabilities as predictors instead of the above text embeddings. Considering there are only four classes in our specification, there is no need to use LASSO. Instead, we use a logistic regression model. However, since the class probabilities are compositional data, we transform the probabilities into additive log ratios, using disruptive conduct as the reference category. We again obtain substantial improvement in out-of-sample prediction accuracy (29.7% for females and 22% for males) when including class probabilities along with the matrix of observed covariates (Figure (ref)).
For each cross-fitting fold, we assess variable importance. For this, McFadden's pseudo $R^{2}$ is computed first for the full model as $1-\frac{\text{deviance all covariates}}{\text{deviance null model}}$. Next we iteratively drop one covariate at a time and compute McFadden's pseudo $R^{2}$ for each model. The reduction in pseudo $R^{2}$ relative to the full model, which includes the class probabilities as predictors, is displayed in Figure (ref). The interpretation of variable importance should account for the reference category (Disruptive Conduct). Here, the reduction is shown as box plots, since we use five-fold cross-validation for prediction, which provides us with a different model for each fold. In addition, panels (c) and (d) provide the marginal effects of increasing community support relative to the disruptive conduct baseline category on predicted probability of recidivism, holding everything else at their median values and fixing the white dummy at white and the education level to its lowest category. We see a sharp reduction in the likelihood of recidivism with increasing community support for both male and female residents relative to disruptive conduct (baseline category)\footnote{The estimates used for this analysis correspond to the model saved for the 5th fold of the 5-fold cross-validation exercise.}. Note by definition, 1 unit increase in the ALR means $e \approx 2.718$ increase in the ratio of community support proportion to disruptive conduct proportion in the classification of messages sent by an individual. Since this is a non-linear model, the changes in the marginal effects for a 1 unit change in ALR vary with the value of ALR and can be represented as the marginal curves in panels (c) and (d). For example, panel (c) shows that as the ratio of community support proportion to disruptive conduct proportion increases from $\exp(-1)\approx 0.36$ to $\exp(0)=1$, the predicted probability of recidivism reduces from around 0.22 to 0.12 for the female unit. Similarly, panel (d) shows that as the ratio of community support proportion to disruptive conduct proportion increases from $\exp(-2)\approx 0.14$ to $\exp(0)=1$, the predicted probability of recidivism reduces from around 0.40 to 0.16 for the male unit.
Further, in Figure (ref)(a), we assess the (out-of-sample) predicted probabilities from our compositional logistic regression model with zero-shot classes against the true recidivism outcome. The comparison of the two histograms shows that the logistic model assigns a higher probability to residents who eventually recidivate. In Figure (ref)(b), we show a calibration plot between the mean predicted probability and the mean observed outcome in 20 quantiles. To do so, we divide the predicted probabilities into 20 quantiles. In each quantile, we compute the mean of the predictions and plot it against the mean of the observed outcomes for individuals in that quantile. If our prediction model is accurate, we expect these points to lie close to a 45-degree line. The figure shows that most points lie close to the red 45-degree line, indicating good calibration of our model's predictions. Similar exercise is done for the male unit, and the results are provided in Figure (ref) in the Appendix.
Finally, we provide prediction results on recidivism aggregated by receiver profiles with and without the class probabilities (Figure (ref)). We find that the prediction accuracy improves by 35% and 23% for the female and male units, respectively.
We compare five estimators for four network models in terms of their ability to replicate observed salient features of the TC networks. These estimators/models are (1) tetra logit estimator (tet.logit) and the (2) joint fixed effect maximum likelihood estimator (jfe) in graham2017econometric that allow for single latent variable (additive) and covariates, 3) Spectral embedding in RDPG model athreya2017statistical which allows for multivariate latent factors but no covariates (rdpg), 4) MLE in latent space model that allows for both multiplicative and additive latent factors but no observed covariates hoff2002latent (lsm) and finally 5) MLE in latent space model from ma2020universal that allows for both multiplicative and additive latent factors and observed covariates (lsm-cov).
Among the above-listed estimators that allow for observed covariates for network formation, we use the time difference (in days) between the residents' entry dates as a node-pair level covariate. The entry dates of residents can be assumed to be random and not explained by their observed or unobserved characteristics, since they are determined by a complex function of court decisions and the availability of beds in TCs. However, the proximity of entry dates is an important determinant of whether residents will exchange messages since overlap in their time in the TC is a necessary condition to exchange messages. We do not use the proximity in exit dates or time in TC as a predictor of network formation, as the unobserved heterogeneity of residents likely impacts these variables, and they are not pre-determined to network formation.
We compare the estimators in terms of their ability to replicate three widely observed structural network properties: modularity girvan2002community, standard deviation of row means, and transitivity or the clustering coefficient newman2018networks. The steps deployed for comparison are as follows. Let us consider the network data for the male unit. We implement each estimator on the observed male unit network to estimate the respective fitted models. Next, we simulate 200 networks using these fitted models (200 for each estimator listed above). The structural properties over these simulated networks are then compared with the true network to assess how well the simulated data can replicate the truth. Figure (ref) displays the performance. We observe that the latent space model that allows for both multiplicative and additive latent factors, along with covariates (proximity in entry dates), outperforms all other models for every property in the male unit. For the female unit, it performs the best for the standard deviation of row means and transitivity of the clustering coefficient, and at least as good as the other estimators for modularity.
This section provides our results on the peer effects model, where the outcome variables are the class probabilities generated by Zero-Shot classification. Our preferred specification involves estimating the network formation model using a latent space model, which allows for both multiplicative and additive latent variables and observed covariates (lsm-cov). This selection is based on the results obtained in the previous section, which show that the latent space model with covariates performs the best in replicating key structural properties of the observed male and female unit networks.
As a first step, we assess the strength of the instruments for explaining the variability in the three endogenous regressors using the Cragg-Donald F statistic (CD statistic) and compare it with the Stock Yogo critical values stock2002testing. Table (ref) shows the CD statistic is between 12 to 14.68 while the critical value at the 10% maximal IV relative bias is 10.25. This indicates the instruments are moderately strong. We note that maximal size-based Stock Yogo critical values are not available when the number of endogenous regressors is greater than 2 (three in our setup, one for each ALR), and hence we compare with only the relative bias critical values.
Tables (ref) and (ref) provide the results for IV 2SLS in text class probabilities after adjusting for the unobserved heterogeneity in the outcome model (Algorithm (ref)) for females and males, respectively. The dimension of the multiplicative factor in the latent space model is chosen using data-driven cross-validation methods for networks in li2020network. The choice of data-driven dimensions for the multiplicative factors is 14 for the female and 19 for the male units, respectively. In addition, the model for network formation contains one additive nodal latent variable and an observed edge-level covariate that captures the difference in the residents' entry dates.
First, we discuss the diagonal elements in Table (ref), which inform us about the spillovers of peer probabilities on female residents for the same class. There is evidence of a peer effect in all three classes relative to disruptive conduct at the 5% significance level. Moreover, there is also evidence of peer spillovers from personal growth and community support categories on rule violations category. We can interpret these results as follows, for example. An individual is more likely to send messages of peer community support if their peers are also sending messages of community support. As a comparison, Table (ref) presents the results when unobserved heterogeneity in the outcome model and network formation are not accounted for.
Next, the same analysis is repeated for the female unit, but the aggregation is done by the receiver profiles (panel b in Table (ref)). When aggregated by the receiver profiles, we find evidence of direct spillover effects only in personal growth ALR and rule violations ALR at the 95% confidence interval. As for the indirect spillovers, we see a positive indirect spillover of peer personal growth on community support. On the contrary, peer community support negatively impacts personal growth. This negative coefficient might emanate from highly positive correlations between peer personal growth and community support ALR relative to disruptive conduct.
The peer effects estimation is next done for the residents in the male units (Table (ref)). The estimates, when aggregated by sender profiles, suggest that there are statistically significant positive direct spillovers of community support and rule violations ALR at the 95$\%$ confidence level. For all of the indirect spillovers, the 95% confidence interval includes 0. In panel (b) of Table (ref), estimates are provided for receiver profile aggregation. The direct spillovers are all positive, but only the ones for personal growth and rule violations are statistically significant at the 5% level. The indirect spillovers are again small relative to the standard errors. For comparison, Table (ref) provides the estimates without adjusting for unobserved heterogeneity.
Overall, we find evidence of peer influence in the ALR of text class probabilities both when aggregated by the sender profiles and when aggregated by the receiver profiles. Interestingly, the magnitude and standard errors change considerably once we account for unobserved confounding in outcome and network formation using latent variables. This underscores the importance of controlling for latent confounding in peer effect studies. We also note from our results that the same category peer influence spillovers generally have a larger impact than across-category spillovers for both male and female units. Finally, we note that even though some of the values in our peer effect parameter matrix are greater than 1, the matrix $(I-D \otimes G)$ is invertible, which makes the structural multivariate peer effect model in Equation (ref) stable (by Lemma 1 in zhu2020multivariate).
\paragraph{Robustness of Peer Effect Estimates:}
We validate the peer effect estimates obtained in Section (ref) by comparing them against the distribution of peer effects obtained by randomly permuting peers. Intuitively, the “significant” peer effects at conventional levels should be much higher in standardized absolute magnitude than this distribution since the peer effects with random peer assignment should be null. This is implemented by randomly shuffling the network connections, preserving each resident's unique number of connections and the total number of messages exchanged. This process is repeated 50 times, yielding 50 simulated random networks. Once these networks are obtained, we repeat the steps outlined in Algorithm (ref) to estimate the peer effects matrix. Note that, since the network changes with each simulation, we re-estimate the latent variables using our preferred latent-space model with covariates. The 50 networks and peer effect estimation can be used to construct the entire empirical distribution of the t-statistic of the peer effects parameter. The final step is to compute the empirical p-value as the proportion of test statistics that exceed the observed statistic value.
The empirical p-values for the observed test statistics of peer effects are reported in Table (ref). The diagonal elements provide same-category peer effects, and the off-diagonal elements provide the cross-category effects. The table shows that each diagonal element has a very small empirical p-value (<0.05). For the male unit, the observed t-statistics for peer effects on Community Support ($p = 0.00$) and Rule Violations ($p = 0.00$) are significant at conventional levels, whereas the peer effect on Personal Growth is not statistically significant ($p = 0.60$). This is consistent with the observations from the main results in Table (ref).
The results in Section (ref) show that there are substantial peer effects in language within the correctional units. The prediction analysis provides evidence that language profiles contain meaningful signals for predicting recidivism, especially in settings where conventional pre-entry covariates have limited predictive power. However, we consider a limitation of our analysis. For policy implications, one has to move beyond predictive power. To establish the value of language profiles for recidivism, it is necessary to estimate the causal link between language profiles and post-release outcomes.
The directed acyclic graph for this causal effect is provided in Figure (ref). The primary objective is to estimate the direct effect of the $i^{th}$ resident's language profile ($L_i$) on the $i's$ recidivism ($R_{i}$), conditional on the observed covariates $X_i$. However, the unobserved heterogeneity $U_{i}$ confounds this relationship. These unobservables can capture multiple factors such as behaviors, personality traits, motivation levels, and aspirations.
We estimate the linear relationship between $R_i$ and $L_i$ controlling for the observed covariates. Next, we do sensitivity analysis using two complementary methods oster2019unobservable and cinelli2020making to get the bounds on how much omitted variable bias is needed relative to the observables to explain away the entire effect of $L_i$ on $R_i$.
First, we implement the approach outlined in oster2019unobservable. oster2019unobservable defines the “coefficient of proportionality" $\delta$ that informs the extent of the selection on unobservable relative to the observed covariates needed to completely explain the estimated effect of $L_i$ on $R_i$. To calculate $\delta$, the researcher has to fix the value of $R_{\text{max}}$, which is the true $R^2$ if we were to regress $R_i$ on $L_i, X_i$, and the unobservable $U_i$. For our analysis we set $R_{\text{max}}=1.3\tilde{R}$ using the recommendations in oster2019unobservable where $\tilde{R}$ is the $R^2$ from the regression of $R_i$ on the observables $L_i,X_i$.
Table (ref) provides the results of this analysis for the female sender, female receiver, male sender, and receiver in panels A--D, respectively. Column (1) reports the estimated $\beta$ from the regression of $R_i$ on $L_i$ controlling for the observed covariates (LSI-R, age, education, and race). In Column (2), we report the $\beta$ for the benchmark scenario where the strength of the selection on unobservables is equivalent to that of the selection on observed covariates. Column (3) reports the $\delta^*$ that will eliminate the effect of $L_i$ on $R_i$. Finally, the last column reports the $R^2$ from the model when $R_i$ is regressed on both the language profiles and the covariates.
For the female unit, all three language profile coefficients remain economically meaningful even when considering that the selection on unobservables is as strong as the selection on observables. For instance, in the female sender profiles, Community Support reduces recidivism by 2.4 percentage points (versus 3.7 in OLS), while Rule Violations increases it by 7.6 percentage points (versus 11.8 in OLS). We see similar patterns, but only for Community support and Rule Violation ALRs in the male unit. Lastly, the results in column (3) suggest that the selection on unobservables has to be in the range of 2.7--2.9 times the selection on observables to completely explain away the impact of language profiles on recidivism.
For the second bounding analysis presented in Table (ref), we use the “sensmakr” method in cinelli2020making. Unlike our previous bounding analysis, sensemakr R package does not support simultaneous sensitivity analysis for multiple correlated treatments. Therefore, we construct a single language treatment variable, which we call the “positive language" variable. To do so, we add the positive zero-shot classes Personal Growth and Community Support to create a positive language class, and the negative zero-shot classes Rule Violations and Disruptive Conduct to create a negative class. We then take the log ratio of these two new classes to obtain our positive language variable.
From the Table (ref), we see the Robustness Values $RV(q=1)$ ranges from 0.22 to 0.24. This indicates that, with the assumption of equal strength of association with the outcome and the treatment (in terms of partial $R^2$, i.e., explaining the residual variance), a confounder needs to explain 22-24% of the residual variation in both the outcome and the treatment in a model that also contains the observed covariates. A similar observation can be made from the column of $RV(\alpha=0.05)$ in terms of the strength of confounders needed to completely invalidate the statistical significance of the effects of positive language. Benchmarking against the observed covariate LSI-R, we see that the statistical significance of the effect of positive language will remain (for female senders) even if we have an unobserved confounder that is 3 times as strong as LSI-R in terms of association with both the outcome and the treatment.
The bounding analysis above provides an estimate of the strength of unobserved confounders required to invalidate our estimate. However, we note that to satisfactorily establish the causal impact of language profiles on recidivism, one needs either to implement a randomized controlled trial that creates random variation in language or to use natural experiments to tease out the causal effect. We see this as a future direction of research, building on our work on the underlying mechanisms of peer effects in language interactions in these facilities.
In conclusion, we develop a framework that combines AI-based measurements of unstructured text with econometric methods to investigate peer effects in language. This framework can be applied more generally, outside our application to recidivism, to study peer effects in contexts such as education and worker productivity, where written or spoken language exchange data are available. Our results also highlight the predictive ability of AI measurements to complement traditional covariates for predicting economic outcomes.
{-1.2pt}