1 Introduction
In many disciplines, including epidemiology, psychometrics, political science, and text mining, longitudinal network models are widely used to study the dynamic interactions among items. Here ”items” have the same meaning as ”nodes” in network models. Examples include friendship dynamics (Knecht, 2008), changes in the relationships among various attributes of hotel text review data (Han et al., 2016), and longitudinal survey data for reporting the experiences of participants in psychology (Geschwind et al., 2011). This has practical implications; for instance, studying how the relationships among the sentiment keywords of hotel reviews change over time can be useful for designing marketing strategies. In this manuscript, we propose a new model to directly interpret the temporal interactions among items for longitudinal network data. We embed timevarying (or functional) interaction parameters within wellestablished exponential random graph model (ERGM) frameworks. Such models face several computational and inferential challenges: (1) The models include intractable normalizing constants because we do not assume conditional independence among the items, and (2) With an increasing number of items, the number of functional parameters grows exponentially. It is challenging to identify significant interaction parameters using standard variable selection approaches because those methods do not consider temporal trends in functional parameters. To address these challenges, we develop a novel Bayes approach based on the functional shrinkage method (Shin et al., 2020)
combined with an auxiliary variable Markov chain Monte Carlo (MCMC) algorithm
(Murray et al., 2006, Liang, 2010). Our method can automatically detect strong temporal interactions among items, while avoiding the direct evaluation of the intractable normalizing constants included in the model.There is an extensive literature on temporal network modeling. Exponential random graph models (ERGMs) (Robins et al., 2007) are widely used to study static networks. By extending ERGMs, Hanneke et al. (2010) develops temporal ERGMs (TERGMs) to analyze networks observed in discrete times. To address parameter interpretation issues in the incidence of ties and their duration in TERGMs, Krivitsky and Handcock (2014) suggest separate temporal ERGMs (STERGMs) that use formation and dissolution networks in dynamic network modeling. The latent space approach (Hoff et al., 2002) is another important class of network models that can measure the similarities between actors using the distance between them in a latent space. Sewell and Chen (2015) develops a dynamic latent space model that allows temporal trajectories for each actor in a latent Euclidean space. Both random graph and latent space type models are widely used to investigate the dynamic relations among nodes from binary network data. In particular, they are useful for studying the overall patterns of the network such as clustering. With somewhat different aspects, there have been a number of recent proposals for network psychometrics to study binary item response data sets (van Borkulo et al., 2014, Jin and Jeon, 2019). These provide an elegant approach to detect the interactions among items from respondents. Such data sets are common, for instance, in binary survey responses (yes, no) for items from individuals, or in binary review data (presence, absence) about aspects from hotels. Here, we focus on studying the temporal network structures among times for such data sets.
Several longitudinal models have also been developed for inferring a dependence structure among items. These are based on latent growth models (cf. Rao, 1958, Tucker, 1958), which represent repeated measures of binary responses as a function of time and other variables under the structural equation modeling framework. Andrade and Tavares (2005) study longitudinal item response data by assuming the conditional independence of the responses to the items. To address the serial dependence in longitudinal categorical item response analysis, Dunson (2003) and Vasdekis et al. (2012) incorporate a firstorder autoregressive structure for the latent variables in a latent growth model, while Segawa (2005)
uses an autoregressive model for the error terms in the structural part of a latent growth model.
Bringmann et al. (2013)uses a multilevel vector autoregressive model that can combine betweensubject and withinsubject information.
Jeon and RabeHesketh (2016) applies an autoregressive structure directly to observed item responses in longitudinal binary item response analysis. Such models can capture the temporal dependence between binary responses to the same items as well as account for heterogeneity in individual responses by adding random effects. However, it is difficult for these models to identify the interactions among items.In this manuscript, we propose functional inhomogeneous exponential random graph models (FIERGMs) to investigate the temporal interactions between items. Our model can detect strong pairwise interactions without making the conditional independence assumption among items. There are several inferential and computational challenges for FIERGMs. The likelihood functions involve intractable normalizing constants, which result in doublyintractable distributions (Murray et al., 2006) in Bayesian analysis. Furthermore, with an increasing number of items (), the number of model parameters grows at order . Since these parameters are timevarying functional objects, the conventional variable selection method can not be directly applicable. Even in static networks, relatively few approaches have been developed for variable selection. Based on the conditional independence assumption among items, van Borkulo et al. (2014) imposes an penalty on Ising graphical models (Ravikumar et al., 2010). However, such an approximation becomes poor in the presence of a strong dependence on the network. To address this, Park and Jin (2019) develops a Bayesian model selection approach for item response data. However, to our knowledge, no existing approaches provide a variable selection procedure for longitudinal networks. This motivates the development of new methods that allow shrinkage for timevarying parameters. To address these challenges, we adopt two methods: (1) a double MetropolisHastings (Liang, 2010)
, an auxiliary variable MCMC method that cancels out the intractable normalizing constants in the acceptance probability, and (2) a functional horseshoe prior
(Shin et al., 2020) that encourages shrinkage toward the zero function for weak signals. We show that our methods can recover the dependence structure of longitudinal networks as well as detect strong timevarying interactions.The rest of this manuscript is organized as follows. In Section 2, we propose functional inhomogeneous exponential random graph models (FIERGMs) and discuss their computational and inferential challenges. In Section 3, we describe the functional horseshoe double MetropolisHastings (FHSDMH) to make an inference for the FIERGMs and describe the implementation details. We provide a guideline to shrink the functional parameters, which can automatically detect the significant temporal interactions. In Section 4, we apply the FHSDMH to three real data sets: Korea youth panel survey data, motivation to succeed survey data, and hotel review data. In Section 5, we investigate the performance of the proposed method through simulation studies. We conclude with a discussion in Section 6.
2 Functional Inhomogeneous Exponential Random Graph Models
Let denote data with responses to binary items observed through time points. For all , if the th individual responses th item positively (or correctly) at time ; otherwise . To account for the pairwise temporal interactions among items, we propose the functional version of IERGMs (Frank and Strauss, 1986) as
(1) 
For time , is an item easiness parameter, which acts as the intercept term for item and is the pairwise interaction among items. Consider the model parameter , where . We can define and as denoting the th row of and th column of , respectively. Then, is the parameter for time and is the th functional parameter over time. Our models assume that (the observed binary response at time ) is only dependent on . In Section 3, we introduce functional priors to account for the temporal trends within each functional parameter .
Inferences for FIERGMs are challenging because of the intractable normalizing constants included in (1). At each time , the calculation of requires summing the overall configurations of the binary response data, which is intractable even with the moderate sizes of and . For such models, frequentist methods (van Borkulo et al., 2014) have been developed based on the pseudolikelihood approximation (Besag, 1974) by assuming conditional independence among items. However, when there is strong dependence in the data, which is often the case in practice, the pseudolikelihood approximation performs poorly.
Another difficulty with FIERGMs is that the number of model parameters increases at an order of ( represents the number of items). To rule out weak interaction parameters, Park and Jin (2019) develops a Bayesian variable selection method for IERGMs in static networks. However, compared with their problem, our case is more complex because of the functional parameters that vary across each time point. To address this challenge, we propose a novel MCMC approach based on a shrinkage prior on function spaces. Our method can automatically detect functional parameters with a weak signal, while providing posterior samples from an intractable likelihood function.
3 Functional Horseshoe Double MetropolisHastings
For FIERGMs, Bayesian approaches are useful to capture the dependence structure in temporal networks without assuming conditional independence among items. Furthermore, we can easily incorporate shrinkage priors to rule out parameters with weak signals. In this section, we propose an MCMC algorithm for FIERGMs. We combine double MetropolisHastings (DMH) (Liang, 2010) with the functional horseshoe prior (Shin et al., 2020) to address the computational and inferential challenges in FIERGMs.
3.1 Bayesian Hierarchical Models with the Functional Horseshoe Prior
Since the number of model parameters () for FIERGMs increases exponentially, one might consider standard shrinkage priors (cf. George and McCulloch, 1993, Carvalho et al., 2010) to rule out weak signals. However, such methods do not account for the temporal trends in the functional parameters in FIERGMs. Furthermore, it is not clear how to select the functional parameters through standard shrinkage priors because they can only induce sparsity on individual . Therefore, in this manuscript, we apply the functional horseshoe (FHS) prior (Shin et al., 2020), which can impose shrinkage on the shape of functions. The FHS can encourage shrinkage toward any parametric class of functions. Here, we focus on shrinkage toward zero functions to detect strong signals, which allows us to perform a natural model selection in FIERGMs. Shin et al. (2020) shows that the posterior constructed from the FHS prior is concentrated at a nearoptimal minmax rate.
Consider the th functional model parameter in (1). Then, the prior on the functional model parameter can be
(2) 
where is a matrix of prespecified basis functions and is a vector of basis coefficients. explains an error that cannot be captured by a mean trend . We assume that the nonparametric basis expansion can capture the temporally dependent trends within each functional model parameter . Here, we use the Bspline basis (De Boor et al., 1978), but other basis functions can also be used. The FHS can shrink toward the null function subspace spanned by a null regressor matrix with
. Then, we can define the FHS hyperpriors as
(3) 
To impose shrinkage toward the zero function, we set the null function space as . Then, and , where . Following Shin et al. (2020), we set and to satisfy the near optimal nonparametric posterior contraction rate. We can summarize the hierarchical models for FIERGMs as
(4) 
where
(5) 
3.2 Markov chain Monte Carlo Implementation
Our model (5) includes intractable normalizing constants , which pose inferential and computational challenges. The resulting posterior (4) is called a doublyintractable distribution having extra unknown normalizing terms , which cannot be canceled out in standard MCMC approaches. Several Bayes methods have been developed for sampling from doubly intractable distributions. Just few of these include constructing Monte Carlo approximations for (cf. Atchade et al., 2008, Lyne et al., 2015, Park and Haran, 2020) and generating auxiliary variables to avoid the direct evaluation of (cf. Murray et al., 2006, Liang, 2010). Considering that constructing Monte Carlo approximations is unstable with an increasing number of parameters, auxiliary variable approaches may be appropriate for the problems considered in this manuscript. In particular, double MetropolisHastings (DMH) (Liang, 2010) is the most practical for computationally expensive problems and the only feasible approach for highdimensional parameter problems among current approaches (see Park and Haran (2018)
for comparisons between algorithms). Therefore, in what follows, we incorporate FHS priors with DMH to impose functional shrinkage as well as estimate the model parameters for FIERGMs.
Consider the model parameter
, and the hyperparameters for the FHS priors
. We update the parameters sequentially through the Gibbs sampler. Let the parameters at the th iteration of the Markov chain be(6) 
For we use the arbitrary initial values with . We can update the parameters successively. The th model parameter for time can be updated from
(7) 
where . Since includes intractable , we use a double MetropolisHastings (DMH) algorithm (Liang, 2010) to update . This is a nested MCMC algorithm; a MetropolisHastings sampler is implemented within another MetropolisHastings sampler. At each iteration of the MCMC (outer MCMC) is proposed from the proposal . For the given , DMH simulates an auxiliary variable from the probability model through the standard MetropolisHastings sampler (inner sampler). For each iteration of the inner sampler, pairs from are randomly chosen; is set to or based on the full conditional probabilities of the networks. See Hunter et al. (2008) for more details. Theoretically, we can simulate an exact auxiliary variable as the inner sampler length approaches infinity; of course the length of the inner sampler should be finite in practice. Following Liang (2010), we use the inner sampler length as , which is equal to the sample size (number of respondents). This generates the reasonably accurate auxiliary variable. Then, the resulting acceptance probability for updating is
(8) 
We note that (8) does not include intractable normalizing constant . The main idea of this approach is to cancel out in the acceptance probability with a clever choice of the auxiliary variable. The more the simulated is close to the observed , the more likely the proposed will be accepted. We repeat this procedure for and .
Then, for can be obtained from
(9) 
We transform to use a slicer sampler for better mixing as well as computational efficiency. The remaining parameters and can be updated from
(10) 
where the conditional distributions are an inverse gamma distribution and a normal distribution, respectively. The conditional distributions for all parameters are described in the supplementary material. The FHSDMH algorithm is summarized in Algorithm
1.
, where .

.

Take .
With each iteration of the MCMC, our algorithm generates an auxiliary variable for times to update (see Step 2 of Algorithm 1). Since we use the inner sampler length as (identical to the number of respondents), the computational complexity of the brute force implementation is , where is the number of items and is the number of observed time points. However, in (5), we assume that is only dependent on . Therefore, we can use parallel computing to generate independently of the given . Then the computational complexity of our method is , where is the number of available processors. The parallel computation is implemented through the OpenMp library in C++. FHSDMH can automatically shrink the parameters having weak signals to zero functions as well as generate posterior samples from the complex hierarchical models in (5).
3.3 Shrinkage Procedure for Functional Parameters in FIERGMs
For FIERGMs, we have the functional parameters , where increases exponentially with the number of items (). Among the functional parameters, the first of them represent the easiness of the corresponding items across time points. These can be regarded as the intercept terms in standard regression models. The remaining number of parameters represent the pairwise interaction among the items across time points. Since not every item has statistically significant interactions, it is important to identify the important interactions only by shrinking the others to zero.
Consider the reparameterization , where can be interpreted as the weight that the posterior mean for function places on the null function subspace (Shin et al., 2020). Therefore, a larger indicates higher weights on the zero function for . According to Shin et al. (2020)[Theorem 3.2], the posterior distribution of converges toward 1 when the true shape of is the zero function. On the contrary, if the true , the posterior distribution of contracts toward 0. Following Shin et al. (2020), we set the threshold for as 0.5. If the posterior mean of is greater than 0.5, we diagnose ; otherwise
4 Applications
Here, we illustrate the application of our method to three real data examples: (1) Korea youth panel survey data, (2) motivation to succeed survey data, and (3) hotel review data. We observe that FHSDMH can shrink weak interactions toward the zero function as well as recover the dependence structure of the longitudinal network well. For the MCMC implementation, we use an independent normal proposal. The convergence of MCMC methods has been checked by the Monte Carlo standard errors
(Jones et al., 2006, Flegal et al., 2008).4.1 Korea Youth Panel Survey Data
Data Description
The first example came from the Korea Youth Panel Survey (Lee et al., 2010) that tracked a nationally representative sample of secondyear middle school students for six consecutive years from 2003 to 2008. Following Jeon and RabeHesketh (2012), we analyzed the subset of the data that excludes about 2% of the students who changed their school membership during their middle school and/or high school years.
In this application, we consider the selfesteem scale that consists of 12 items on a 5point Likerttype scale; the response options are “Strongly disagree”, “Disagree”, “Neither agree nor disagree”, “Agree”, and “Strongly agree”. We choose the seven items analyzed by Jeon and RabeHesketh (2016): (1) I sometimes think I am a useless person, (2) I sometimes think I am a bad person, (3) I sometimes feel like I am a failure, (4) I think I am a trouble maker, (5) I think I am a juvenile delinquent, (6) Other people think I am a trouble maker, and (7) Other people think I am a juvenile delinquent. To measure a positive selfimage (or selfesteem), the response categories were reversed and dichotomized at point 3 ( 3 and 3). The proportion of male students was 49.9%. At each time point, the data include binary responses from about individuals to the items; this results in functional parameters in FIERGMs.
Interactions  Time 1  Time 2  Time 3  Time 4  Time 5  Time 6 

0.756  0.852  0.767  0.889  0.899  0.751  
0.648  0.473  0.631  0.707  0.609  0.775  
0.597  0.447  0.517  0.351  0.616  0.435  
0.705  0.756  0.782  0.752  0.549  1.011  
0.400  0.450  0.555  0.584  0.494  0.353  
0.667  0.524  0.501  0.344  0.566  0.672  
0.908  0.822  0.715  1.033  0.828  0.825  
0.678  0.729  0.329  0.775  0.527  0.553  
3.050  3.048  3.192  2.890  3.319  3.352  
0.842  0.758  0.684  0.803  0.885  0.951  
0.805  0.585  0.800  0.593  0.891  0.779  
0.558  0.365  0.446  0.524  0.388  0.547 
Analysis Results
Among the 28 functional parameters, FHSDMH shrinks 10 of them to zero functions (Figure 2 (B)). Figure 1 and Table 1 describe the estimated dependence structures and their nonzero interaction parameters. We observe that every connection of items (nodes) shows positive relationships and that the item dependence structure is consistent over time. This indicates that these items measure a construct in the same direction and that the structure of the construct (or the relationships between the construct and items) is stable over 6 years.
Items (4) and (6) have the strongest positive interaction, which makes sense given that both items are about images as trouble makers. Interestingly, item (4) is connected to all the items except item (7). This may indicate that having a positive selfimage as not a trouble maker is important to have an overall positive view of oneself, while it does not seem to strongly impact how one thinks about others’ view of oneself as a juvenile delinquent.
Two additional patterns are identified in the estimated temporal network (Figure 1). The items related to selfevaluation (items (1)(5)) are connected to one another. Similarly, items (6) and (7), which are related to reputation from others, are connected. Among the selfevaluation items, only item (5) is connected to items (6) and (7) representing reputation from others. This connection indicates that students who regard themselves as not juvenile delinquents tend to be more sensitive to others’ views of themselves.
Model Validation
To validate our method, we compare the summary statistics between the observed network and fitted network from the posterior predictive distribution
(Gelman et al., 2013). For observed data , we use the following summary statisticsIn our example, because we have 28 parameters for 6 time points. For a given posterior sample from FHSDMH, we simulate binary response data . Then we obtain summary statistics . If these synthetic summary statistics resemble the observed summary statistics well, then our FHSDMH posterior sample can be regarded as a reasonable approximation of the true posterior distribution. To implement this, we obtain 1,000 thinned posterior samples from FHSDMH. We simulate from the 1,000 posterior samples. Then, we calculate the sample mean of the summary statistics . In Figure 2 (A), the mean of the observed statistics and simulated statistics follow a straight line, which indicates that our fitted model is a good reflection of the true network structure of the observed data.
4.2 Motivation to Succeed Data
Data Description
The data for the second application came from the Pathways to Desistance study (Mulvey et al., 2004)
, which is a multisite longitudinal study that follows 1,354 serious juvenile offenders from adolescence to young adulthood in Philadelphia and Phoenix. Participants completed baseline interviews in November 2000 and followup interviews at 6, 12, 18, 24, 30, 36, 48, 60, 72, and 84 months postbaseline (first follow up interview completed in May 2001; last follow up interview completed in March 2010). The aims of the study are to identify initial patterns of how serious adolescent offenders stop antisocial activity; describe the role of the social context and developmental changes in promoting these positive changes; and compare the effects of sanctions and interventions in promoting these changes.
From this large study, we used the motivation to succeed scale (Eccles et al., 1998), which includes six items that measure the respondent’s assessment of the opportunities available in his/her neighborhood regarding schooling and work. An additional two items were included that measure the adolescent’s perceptions of how far he/she would like to go in school and how far he/she think they will go in school. The eight test items are as follows: (1) In my neighborhood, it is easy for a young person to get a good job; (2) Most of my friends will graduate from high school; (3) In my neighborhood, it is hard to make money without doing something illegal; (4) College is too expensive for most people in my neighborhood; (5) We have fewer opportunities to succeed than kids from other neighborhoods; (6) Our chances of getting ahead/being successful are not very good; (7) How far would you like to go in school? (8) How far do you think you will go in school?
The responses to the first six questions are on a fivepoint Likertscale (“Strongly disagree”, “Disagree”, “Neither agree nor disagree”, “Agree”, and “Strongly agree”). The responses to the last two questions are “Drop out before graduation”, “Graduate from HS”, “Go to business, tech school or jr college”, “Graduate from college”, and “Go to graduate or professional school”. We dichotomized the responses to the first six items by assigning 1 to “Agree”, and “Strongly Agree” responses and 0 otherwise. For the last two items, we assigned 1 to “Graduate from college” and “Go to graduate or professional school” and 0 otherwise. After removing cases with missing responses, 740 respondents remained. In summary, at each time point, the data include respondents for the items, which results in 36 functional parameters in FIERGMs.
Interactions  Time 1  Time 2  Time 3  Time 4  Time 5  Time 6 

0.444  0.632  0.988  0.858  0.963  1.083  
0.893  0.998  0.878  1.136  1.112  1.263  
1.399  1.229  1.724  1.735  1.643  1.711  
0.914  1.339  1.505  1.221  1.304  1.362  
0.426  0.226  0.620  0.898  0.632  0.486  
2.335  2.597  2.862  2.163  2.388  2.913  
1.663  0.978  0.935  1.529  1.567  0.514  
4.398  4.947  5.465  5.717  5.138  5.997  
Interactions  Time 7  Time 8  Time 9  Time 10  Time 11  
0.982  0.898  0.910  0.505  0.720  
0.783  1.232  1.414  1.773  0.979  
1.836  2.531  2.989  3.008  3.211  
1.371  1.465  0.310  1.045  1.663  
1.278  2.099  4.202  4.545  4.126  
2.928  2.286  2.963  3.436  2.030  
2.027  3.764  3.654  7.281  7.031  
5.945  6.674  7.710  8.440  8.256 
Analysis Results and Model Validation
Figure 4 (B) shows that among the 36 functional parameters, FHSDMH shrinks 22 of them to zero functions. Figure 3 and Table 2 show the estimated network structures and their estimated nonzero interaction parameters. We observe three negative connections and five positive connections among the items. Such relationships are generally consistent over time, whereas their strengths vary to some degree. In particular, item (1) and item (3) are negatively connected. This indicates that students who believe that it is easy for young people to get a job in their neighborhood tend to disagree that doing illegal things is necessary to earn money. On the contrary, item (3) shows a positive relationship with item (4). Students who deem it necessary to engage in illegal activities to earn money also think that college tuition is too expensive. From these findings, we can conclude that poverty is an important factor behind students’ antisocial activities.
In addition, item (6), which represents the chance of success, has negative relationships with items (8) and (2), but a positive relationship with item (5). Such connections show that students who anticipate them having little chance of success are less likely to think that they will opt for higher education. Instead, they think their neighborhoods, including themselves, have fewer opportunities to be successful. These relationships convey that students who think they have insufficient opportunities than others are less likely to think that they will go college or graduate high school; rather, they think their neighborhood has fewer opportunities for them to be successful.
In particular, the strength of the interaction between item (6) and item (2) becomes stronger over time. It can be inferred from this result that earning a high school diploma or higher degree brings students closer to success. Therefore, education is an important factor in leading students to success. Additionally, items (7) and (8) have the strongest relationship; this indicates that a desire and willingness to go to college are closely related. This suggests that the desire to go to college motivates students to overcome the obstacles presented by their living situation, should those obstacles exist, and pursue their dream of higher education.
As in the previous example, we simulate the summary statistics from 1,000 thinned posterior samples obtained from FHSDMH. Figure 4 (A) shows that the mean of the 1,000 simulated summary statistics aligns with the observed summary statistics, indicating that our model fits the observed data well.
4.3 Hotel Review Data
Data Description
Item  Aspectbased Sentiments  Sentiment Keywords 

1  Price  satisfaction  reasonable price, inexpensive price 
2  Price  dissatisfaction  very expensive, terrible price, ridiculous price 
3  Room  satisfaction  amazing room, luxury room, quiet room 
4  Room  dissatisfaction  dry room, dirty room, old room 
5  Subsidiary facilities  satisfaction  clean lounge, comfortable lounge 
6  Subsidiary facilities  dissatisfaction  dirty floor, old lounge 
7  Food  satisfaction  delicious buffet, nearby restaurant, clean cafeteria 
8  Food  dissatisfaction  dirty restaurant, tasteless cafeteria 
9  Interior design  satisfaction  nice place, luxurious interior 
10  Interior design  dissatisfaction  poor lighting, no sunlight 
11  Service  satisfaction  great service, prompt service 
12  Service  dissatisfaction  abysmal service, satisfactory service, surly service 
13  Bed  satisfaction  get a good night’s sleep 
14  Bed  dissatisfaction  terrible bed, uncomfortable bed 
Many reviewers express their opinions by writing reviews on websites, and their satisfaction toward hotels is summarized with sentiment keywords. Sentiment keywords represent emotional expressions in hotel review data; the positive and negative opinions in the review data correspond to specific keywords. We use text mining to construct aspectbased sentiments, which are composed of keywords with similar aspects. Here, we study the temporal networks among aspectbased sentiments about hotel reviews from reservation websites. We collect 231,862 reviews of 423 hotels in South Korea from 2018 to 2019 and exclude advertisements and spam reviews. To obtain the sentiment keywords from the data, we use the natural language process through the Daumsoft Text Mining Engine Version 2. With these keywords, we construct 14 aspectbased sentiments such as prices, rooms, subsidiary facilities, food, interior design, service, and bed. Table
3 summarizes the aspectsbased sentiments and their keywords. The data set contains binary review information from hotels over the seven time points from the first quarter of 2018 to the third quarter of 2019. The binary value indicates whether an individual hotel has a review containing aspectbased sentiments (1 existence of sentiments in the review and 0 otherwise). The aspectbased sentiments result in FIERGMs with 105 functional parameters.Interactions  Time 1  Time 2  Time 3  Time 4  Time 5  Time 6  Time 7 

1.651  1.751  1.613  2.142  1.812  2.017  1.809  
1.277  1.317  1.259  1.146  1.818  1.501  1.108  
1.387  0.771  0.727  1.290  1.369  2.039  1.674  
1.411  1.224  1.753  1.489  1.470  1.027  0.861  
1.229  1.176  1.293  1.510  1.821  1.126  0.963  
1.400  0.989  0.890  1.162  1.211  1.634  1.698  
2.222  1.324  1.380  0.885  1.055  0.900  1.166  
0.709  1.463  1.316  1.392  1.278  0.996  1.181  
1.295  1.120  1.198  0.985  1.326  0.879  1.079  
0.766  0.962  0.841  0.698  0.609  0.832  2.220 
Analysis Results and Model Validation
Analysis of Dyadic Relationships
Figure 6 (B) shows that among the 105 functional parameters, FHSDMH shrinks 45 of them to zero functions. Figure 5 and Table 4 describe the estimated network structures and their nonzero interaction parameters. In Figure 5, we changed the layout to the method of Fruchterman and Reingold (Fruchterman and Reingold, 1991), which exploits analogies between the relational structure in graphs, to visualize the node cluster structure efficiently. We provide a circular layout for the hotel review data in the supplementary material. We observe several meaningful patterns. First, item (8) (food dissatisfaction) and item (9) (interior design satisfaction) are connected positively at each time point. This indicates that although reviewers may be satisfied with the interior and exterior design of the hotel (including the restaurant), they are not necessarily satisfied with the taste of the food. On the contrary, satisfaction with the food (item (7)) is linked to satisfaction with the service (item (11)). Since this also has a strong positive relationship at all times, it can be seen as related to satisfaction with the food and service.
At Time 1, there is one negative relationship between room satisfaction (item (3)) and bed dissatisfaction (item (14)). We can thus infer that there are some favorable reviews about room condition but discontent with the bed quality for reasons such as a hard bed, dirty bed, etc. However, the items connecting favorable and unfavorable sentiments tend to display positive relationships because of the characteristics of the review data (i.e., reviewers mention negative and positive sentiment keywords on different aspects simultaneously in reviews). Therefore, the input matrix of the review data according to favorable and unfavorable aspectbased sentiments is analyzed with cooccurrence information. When this occurs, our model labels the interaction between satisfaction and dissatisfaction as positive. In particular, Table 4 shows that food satisfaction (item (7)) and food dissatisfaction (item (8)) have a strong positive relationship at all times, even though they have the opposite meanings. From this, we can infer that reviewers tended to say both good things and bad things about their food.
Analysis of Triangle Relationships
There is a consistent triangle (cyclic) relationship between favorable aspectbased items over time. At Time 1, the connection between room satisfaction (item (3)) and interior design satisfaction (item (9)) is the strongest and their interaction parameter is the highest among the seven time points. Furthermore, both those items are connected to satisfaction with the service (item (11)); hence the three items have a triangle relationship that is consistent at each time point. In other words, many reviews contain a favorable impression of the room, interior design, and service. Furthermore, there is a strong positive relationship between satisfaction with the food (item (7)) and satisfaction with the service (item (11)) at Time 3. In addition, these two items display a strong positive relationship with room satisfaction (item (3)) at all times. Given that items (3) and (11) are connected to item (5) (subsidiary facilities satisfaction), they form another triangle relationship. We can combine this result to conclude that hotel users have favorable opinions about their room, the service, the interior design, the food, and the subsidiary facilities. Considering that such triangle relationships are composed of items (3) and (11), we can infer that satisfaction with the service and with the room conditions are major factors in getting favorable reviews.
As shown in Figure 5 and Table 4, the triangle relationships between items do not change structurally between time points; rather, the strength of these relations vary over time. In terms of unfavorable aspectbased factors, price dissatisfaction (item (1)), service dissatisfaction (item (12)), and meal dissatisfaction (item (8)) are all connected. Further, room dissatisfaction (item (4)), interior design dissatisfaction (item (10)), and bed dissatisfaction (item (14)) are connected. In addition, as shown by the connection strength, the main complaint is that the room design is so poor that it reflects negatively on reviews about the room and bed. Furthermore, the connection between room satisfaction (item (3)), interior design satisfaction (item (9)), and bed satisfaction (item (13)) is maintained over time. Therefore, we can conclude that the interior design influences the overall impression of the room and bed. Other factors are affected by the price such as dissatisfaction with the room and facilities. That is, there is constant connectivity between price dissatisfaction (item (1)), room dissatisfaction (item (4)), and facilities dissatisfaction (item (6)). As in the previous examples, Figure 6 (A) shows that the observed summary statistics and simulated summary statistics are aligned (i.e., our model fits well).
5 Simulated Data Example
To validate our method, we conduct a simulation study under different scenarios. We simulate temporal binary response data sets with observations, items, and time points. This results in 55 functional parameters. We provide cyclic trends for the parameters using a function , which has fluctuating patterns around 0. We simulate 10 intercept parameters () from . Among the 45 interaction parameters (), 22 of them are simulated from ; these are assumed to be the true zero functions. Then, 11 of the interaction parameters are generated from and the remaining 12 parameters are generated from . Here, we consider four noise strength settings (). For the given model parameters , we simulate the data sets via iterations of MetropolisHastings updates (Hunter et al., 2008). Figure 7 illustrates the simulated functional parameters under the different scenarios.
To study the performance of our method, we calculate the true positive rate (diagnose zero for the true zero functions) and true negative rate (diagnose nonzero for the true nonzero functions) for each scenario. Furthermore, we define the root mean square error (RMSE) as
where is the functional parameter estimate obtained from FHSDMH and is the true parameter. Table 5 summarizes the results. Based on the RMSE, we observe that FHSDMH estimates are reasonably close to the true functional parameters. Furthermore, FHSDMH can detect the true zero function well for reasonable noise settings (). With increasing noise () in the simulated data, it becomes difficult to detect the true zero functions because the noises overwhelm the signals from each function.
Scenario  RMSE  TP  TN 

0.11  0.91  0.91  
0.13  0.86  0.88  
0.21  0.68  0.70  
0.27  0.68  0.55 
6 Discussion
In this manuscript, we embed functional parameters in inhomogeneous exponential random graph models to study the temporal interactions among the items. Our models include intractable normalizing constants, and the number of functional parameters increases with an increasing number of items. We combine a double MetropolisHastings algorithm and a functional shrinkage method to address these computational and inferential challenges. Our study to real and simulated data examples shows that FHSDMH can rule out weak interactions among items as well as provide a direct interpretation of temporal trends. Without depending on a conditional independence assumption, our method can recover the dependence structure of the longitudinal networks well. To our knowledge, this is the first attempt to use functional shrinkage for models with intractable normalizing constants, which is an important contribution.
We model the dynamic trends in the functional parameters using the nonparametric basis expansion in (2); the temporal dependencies for the network parameters are specified in the mean (firstorder) function. One might want to model the dynamic structure through the covariance (secondorder) function, such as the AR(1) correlation between and . However, it is practically and computationally challenging to deploy a dynamic secondorder structure on functional shrinkage. For this reason, our methods use a functional horseshoe prior (Shin et al., 2020) by specifying the firstorder mean structure. We expect the nonparametric basis expansion to be capable of detecting temporally dependent trends. Similar to the variable selection methods (van Borkulo et al., 2014, Park and Jin, 2019) for static networks, our methods are suited to data sets with a sufficient number of respondents () and a moderate number of items (), as examples illustrated in this manuscript. Otherwise, it will suffer from a small , large
issues, which can lead to an unreliable inference. As a simple heuristic, we recommend applying our methods to problems with
.The computational methods developed here allow researchers in many disciplines to study temporal interactions among items for binary response data sets. Our methods could be applicable to a variable selection in existing temporal network models as well as for a broader class of functional models. Examples include temporal exponential random graph models (Hanneke et al., 2010) and their variants (Krivitsky and Handcock, 2014), and functional regression models (Ramsay and Silverman, 2007).
Acknowledgement
Jaewoo Park was partially supported by the Yonsei University Research Fund of 2019220194 and the National Research Foundation of Korea (NRF2020R1C1C1A01003868). Ick Hoon Jin was partially supported by the Yonsei University Research Fund of 2019220210 and the National Research Foundation of Korea (NRF2020R1A2C1A01009881). The authors are grateful to anonymous reviewers for their careful reading and valuable comments.
Supplementary Material for Bayesian Shrinkage for Functional Network Models with Intractable Normalizing Constants
Jaewoo Park, Yeseul Jeon, Minsuk Shin, Minjeong Jeon, and Ick Hoon Jin
Appendix A Full Conditionals for Markov Chain Monte Carlo
The conditional distribution for is
(11) 
Let . Using variable transformation
(12) 
Then, we can sample from using the slice sampler as

, where .

. Take .
The conditional distribution of can be derived as
(13) 
Therefore, the conditional distribution of follows an inverse gamma with parameters .
The conditional distribution of is
(14) 
Since in our setting, this follows a normal distribution with mean
and variance
.Appendix B Circular Layout for the Estimated Networks for a Hotel Keyword Data Set
References

Andrade and Tavares (2005)
Andrade, D. F. and Tavares, H. R. (2005).
Item response theory for longitudinal data: population parameter
estimation.
Journal of multivariate analysis
, 95(1):1–22.  Atchade et al. (2008) Atchade, Y., Lartillot, N., and Robert, C. P. (2008). Bayesian computation for statistical models with intractable normalizing constants. Brazilian Journal of Probability and Statistics, 27:416–436.
 Besag (1974) Besag, J. (1974). Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. Series B (Methodological), 36:192–236.
 Bringmann et al. (2013) Bringmann, L. F., Vissers, N., Wichers, M., Geschwind, N., Kuppens, P., Peeters, F., Borsboom, D., and Tuerlinckx, F. (2013). A network approach to psychopathology: new insights into clinical longitudinal data. PloS one, 8(4):e60188.
 Carvalho et al. (2010) Carvalho, C. M., Polson, N. G., and Scott, J. G. (2010). The horseshoe estimator for sparse signals. Biometrika, 97(2):465–480.
 De Boor et al. (1978) De Boor, C., De Boor, C., Mathématicien, E.U., De Boor, C., and De Boor, C. (1978). A practical guide to splines, volume 27. springerverlag New York.
 Dunson (2003) Dunson, D. B. (2003). Dynamic latent trait models for multidimensional longitudinal data. Journal of the American Statistical Association, 98(463):555–563.
 Eccles et al. (1998) Eccles, J., Wigfield, A., and Schiefele, U. (1998). Motivation in: Eiesenberg, n.(ed.), hand book of child psychology, vol. 3, pp: 10171095.
 Flegal et al. (2008) Flegal, J. M., Haran, M., and Jones, G. L. (2008). Markov chain Monte carlo: Can we trust the third significant figure? Statistical Science, 23:250–260.
 Frank and Strauss (1986) Frank, O. and Strauss, D. (1986). Markov graphs. Journal of the American Statistical Association, 81(395):832–842.
 Fruchterman and Reingold (1991) Fruchterman, T. M. J. and Reingold, E. M. (1991). Graph drawing by force‐directed placement. Software: Practice and Experience, 21(11):1129–1164.
 Gelman et al. (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013). Bayesian data analysis. Chapman and Hall/CRC.
 George and McCulloch (1993) George, E. and McCulloch, R. (1993). Variable selection via gibbs sampling. Journal of the American Statistical Association, 88:881–889.
 Geschwind et al. (2011) Geschwind, N., Peeters, F., Drukker, M., van Os, J., and Wichers, M. (2011). Mindfulness training increases momentary positive emotions and reward experience in adults vulnerable to depression: a randomized controlled trial. Journal of consulting and clinical psychology, 79(5):618.
 Han et al. (2016) Han, H. J., Mankad, S., Gavirneni, N., Verma, R., et al. (2016). What guests really think of your hotel: Text analytics of online customer reviews.
 Hanneke et al. (2010) Hanneke, S., Fu, W., Xing, E. P., et al. (2010). Discrete temporal models of social networks. Electronic Journal of Statistics, 4:585–605.
 Hoff et al. (2002) Hoff, P. D., Raftery, A. E., and Handcock, M. S. (2002). Latent space approaches to social network analysis. Journal of the american Statistical association, 97(460):1090–1098.
 Hunter et al. (2008) Hunter, D. R., Handcock, M. S., Butts, C. T., Goodreau, S. M., and Morris, M. (2008). ergm: A package to fit, simulate and diagnose exponentialfamily models for networks. Journal of statistical software, 24(3).

Jeon and RabeHesketh (2012)
Jeon, M. and RabeHesketh, S. (2012).
Profilelikelihood approach for estimating generalized linear mixed models with factor structures.
Journal of Educational and Behavioral Statistics, 37(4):518–542.  Jeon and RabeHesketh (2016) Jeon, M. and RabeHesketh, S. (2016). An autoregressive growth model for longitudinal item analysis. psychometrika, 81(3):830–850.
 Jin and Jeon (2019) Jin, I. H. and Jeon, M. (2019). A doubly latent space joint model for local item and person dependence in the analysis of item response data. Psychometrika, 84(1):236–260.
 Jones et al. (2006) Jones, G. L., Haran, M., Caffo, B. S., and Neath, R. (2006). Fixedwidth output analysis for Markov chain Monte Carlo. Journal of the American Statistical Association, 101(476):1537–1547.
 Knecht (2008) Knecht, A. B. (2008). Friendship selection and friends’ influence. Dynamics of networks and actor attributes in early adolescence. Utrecht University.
 Krivitsky and Handcock (2014) Krivitsky, P. N. and Handcock, M. S. (2014). A separable model for dynamic networks. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):29–46.
 Lee et al. (2010) Lee, K., Lim, H., and Ahn, S. (2010). Korea youth panel study.
 Liang (2010) Liang, F. (2010). A double Metropolis–Hastings sampler for spatial models with intractable normalizing constants. Journal of Statistical Computation and Simulation, 80(9):1007–1022.

Lyne et al. (2015)
Lyne, A.M., Girolami, M., Atchade, Y., Strathmann, H., and Simpson, D. (2015).
On Russian roulette estimates for Bayesian inference with doublyintractable likelihoods.
Statistical science, 30(4):443–467.  Mulvey et al. (2004) Mulvey, E. P., Steinberg, L., Fagan, J., Cauffman, E., Piquero, A. R., Chassin, L., Knight, G. P., Brame, R., Schubert, C. A., Hecker, T., et al. (2004). Theory and research on desistance from antisocial activity among serious adolescent offenders. Youth Violence and Juvenile Justice, 2(3):213–236.

Murray et al. (2006)
Murray, I., Ghahramani, Z., and MacKay, D. J. C. (2006).
MCMC for doublyintractable distributions.
In
Proceedings of the 22nd Annual Conference on Uncertainty in Artificial Intelligence (UAI06)
, pages 359–366. AUAI Press.  Park and Haran (2018) Park, J. and Haran, M. (2018). Bayesian inference in the presence of intractable normalizing functions. Journal of the American Statistical Association, 113(523):1372–1390.
 Park and Haran (2020) Park, J. and Haran, M. (2020). A function emulation approach for doubly intractable distributions. Journal of Computational and Graphical Statistics, 29(1):66–77.
 Park and Jin (2019) Park, J. and Jin, I. H. (2019). Bayesian model selection for ultrahighdimensional doublyintractable distributions with an application to network psychometrics. arXiv preprint arXiv:1911.07142.
 Ramsay and Silverman (2007) Ramsay, J. O. and Silverman, B. W. (2007). Applied functional data analysis: methods and case studies. Springer.
 Rao (1958) Rao, C. R. (1958). Some statistical methods for comparison of growth curves. Biometrics, 14(1):1–17.

Ravikumar et al. (2010)
Ravikumar, P., Wainwright, M. J., and Lafferty, J. (2010).
Highdimensional Ising model selection using
regularized logistic regression.
The Annals of Statistics, 38:1287–1319.  Robins et al. (2007) Robins, G., Pattison, P., Kalish, Y., and Lusher, D. (2007). An introduction to exponential random graph (p*) models for social networks. Social networks, 29(2):173–191.
 Segawa (2005) Segawa, E. (2005). A growth model for multilevel ordinal data. Journal of Educational and Behavioral Statistics, 30(4):369–396.
 Sewell and Chen (2015) Sewell, D. K. and Chen, Y. (2015). Latent space models for dynamic networks. Journal of the American Statistical Association, 110(512):1646–1657.
 Shin et al. (2020) Shin, M., Bhattacharya, A., and Johnson, V. E. (2020). Functional horseshoe priors for subspace shrinkage. Journal of the American Statistical Association. (in press).
 Tucker (1958) Tucker, L. R. (1958). Determination of parameters of a functional relation by factor analysis. Psychometrika, 23(1):19–23.
 van Borkulo et al. (2014) van Borkulo, C. D., Borsboom, D., Epskamp, S., Blanken, T. F., Boschloo, L., Schoevers, R. A., and Waldorp, L. J. (2014). A new method for constructing networks from binary data. Scientific Reports, 4.
 Vasdekis et al. (2012) Vasdekis, V. G., Cagnone, S., and Moustaki, I. (2012). A composite likelihood inference in latent variable models for ordinal longitudinal responses. Psychometrika, 77(3):425–441.
Comments
There are no comments yet.