Missing data

Standard implementations of wavelet trasforms require time series consisting of measurements taken at evenly spaced times, with no missing data. Most functions provided in wsyn make these same assumptions and throw an error if data are missing. The user is left to decide on and implement a reasonable way of filling missing data. Measures of synchrony can be influenced by data filling techniques that are based on spatial interpolation. We therefore recommend that spatially informed filling procedures not be used. We have previously used the simple approach of replacing missing values in a time series by the median of the non-missing values in the time series (Sheppard et al. 2016). This approach, and other related simple procedures (Sheppard et al. 2016), seem unlikely to artefactually produce significant synchrony, or coherence relationships with other variables, but rely on the percentage of missing data being fairly low and may obscure detection of synchrony or significant coherence relationships if too many data are missing. For applications which differ meaningfully from the prior work for which the tools of this package were developed (e.g., Sheppard et al. (2016); Sheppard, Reid, and Reuman (2017); Sheppard et al. (2019); Walter et al. (2017); Anderson et al. (2018)), different ways of dealing with missing data may be more appropriate.

De-meaning, detrending, standardizing variance, and normalizing

A function cleandat is provided that performs a variety of combinations of data cleaning typically necessary for analyses implemented in wsyn, including de-meaning, linear detrending, standardization of time series variance, and Box-Cox transformations to normalize marginal distributions of time series. Most functions in wsyn assume at least that means have been removed from time series, and throw an error if this is not the case. Approaches based on Fourier surrogate (section ) require time series with approximately normal marginals.

The wavelet phasor mean field

The wavelet phasor mean field (the wpmf function in wsyn) depicts the time and timescale dependence of phase synchrony of a collection of time series. If x_n(t), n = 1, …, N, t = 1, …, T are time series of the same variable measured in N locations at the same times, and if W_n, σ(t) is the wavelet transform of x_n(t) and $w_{n,\sigma}(t)=\frac{W_{n,\sigma}(t)}{|W_{n,\sigma}(t)|}$ has only the information about the complex phases of the transform (unit-magnitude complex numbers such as these are called ), then the wavelet phasor mean field is For combinations of t and σ for which oscillations at time t and timescale σ in the time series x_n(t) have the same phase (they are phase synchronized), the phasors w_n, σ(t) will all point in similar directions in the complex plane, and their sum will be a large-magnitude complex number. For combinations of t and σ for which oscillations at time t and timescale σ in the time series x_n(t) have unrelated phases (they are not phase synchronized), the phasors w_n, σ(t) will all point in random, unrelated directions in the complex plane, and their sum will be a small-magnitude complex number. Therefore plotting the magnitude of the wavelet phasor mean field against time and timescale quantifies the time and timescale dependence of phase synchrony in the x_n(t). The wavelet phasor mean field, being a mean of phasors, always has magnitude between 0 and 1.

We provide an example based on supplementary figure 1 of Sheppard et al. (2016). A related technique, the wavelet mean field (section ), was used there, but the wavelet phasor mean field also applies and is demonstrated here. We construct data consisting of time series measured for 100 time steps in each of 11 locations. The time series have three components. The first component is a sine wave of amplitude 1 and period 10 years for the first half of the time series, and is a sine wave of amplitude 1 and period 5 for the second half of the time series. This same signal is present in all 11 time series and creates the synchrony among them.

times1<-0:50
times2<-51:100
times<-c(times1,times2)
ts1<-c(sin(2*pi*times1/10),sin(2*pi*times2/5))+1.1

The second component is a sine wave of amplitude 1 and period 3 years that is randomly and independently phase shifted in each of the 11 time series. The third component is white noise, independently generated for each time series.

dat<-matrix(NA,11,length(times))
for (counter in 1:dim(dat)[1])
{
  ts2<-3*sin(2*pi*times/3+2*pi*runif(1))+3.1
  ts3<-rnorm(length(times),0,1.5)
  dat[counter,]<-ts1+ts2+ts3    
}
dat<-cleandat(dat,times,1)$cdat

The second and third components do not generate synchrony, and obscure the synchrony of the first component. As a result, synchrony cannot be readily detected by visually examining the time series:

plot(times,dat[1,]/10+1,type='l',xlab="Time",ylab="Time series index",ylim=c(0,12))
for (counter in 2:dim(dat)[1])
{
  lines(times,dat[counter,]/10+counter)
}

Nor can synchrony be readily detected by examining the 55 pairwise correlation coefficients between the time series, which are widely distributed and include many values above and below 0:

cmat<-cor(t(dat))
diag(cmat)<-NA
cmat<-as.vector(cmat)
cmat<-cmat[!is.na(cmat)]
hist(cmat,30,xlab="Pearson correlation",ylab="Count")

But the wavelet phasor mean field sensitively reveals the synchrony and its time and timescale structure:

res<-wpmf(dat,times,sigmethod="quick")
plotmag(res)

The wpmf function implements assessment of the statistical significance of phase synchrony in three ways, one of which is demonstrated by the contour lines on the above plot (which give a 95% confidence level by default - the level can be changed with the sigthresh argument to the plotmag method for the wpmf class). The method of significance testing the wavelet phasor mean field plot is controlled with the sigmethod argument to wpmf, which can be quick (the default), fft or aaft. The quick method compares the mean field magnitude value for each time and timescale separately to a distribution of magnitudes of sums of N random, independent phasors.

Each time/timescale pair is compared independently to the distribution, and the multiple testing problem is not accounted for, so some time/timescales pairs will come out as showing “significant” phase synchrony by chance, i.e., false-positive detections of phase synchrony can occur. For instance, the small islands of significant synchrony at timescale approximately 5 and times about 15 and 40 on the above plot are false positives.

Signficance is based on stochastic generation of magnitudes of sums of random phasors, so significance contours will differ slightly on repeat runs. Increasing the number of randomizations (argument nrand to wpmf) reduces this variation.

The quick method can be inaccurate for very short timescales. The two alternative methods, fft and aaft, mitigate this problem but are substantially slower. The fft and aaft methods are based on surrogate datasets (section ) so are discussed in section .

The power and plotphase methods also work on wpmf objects.

Examples of the wavelet phasor mean field technique applied to real datasets are in, for instance, Sheppard et al. (2019) (their Fig. S1) and Anderson et al. (2018) (their Fig. 4).

The wavelet mean field

The wavelet mean field (the wmf function in wsyn) depicts the time and timescale dependence of synchrony of a collection of time series x_n(t) for n = 1, …, N and t = 1, …, T, taking into account both phase synchrony and associations through time of magnitudes of oscillations in different time series at a given timescale. See Sheppard et al. (2016) for a precise mathematical definition. The plot is similar in format to a wavelet phasor mean field plot, but without significance contours:

res<-wmf(dat,times)
plotmag(res)

The wavelet mean field is a more useful technique than the wavelet phasor mean field insofar as it accounts for associations of magnitudes of oscillation, in addition to phase synchrony, but it is less useful insofar as significance contours are not available. The wavelet mean field also has some mathematical advantages, described in Sheppard et al. (2016) and Sheppard et al. (2019). The “wavelet Moran theorem” and other theorems described in those references use the wavelet mean field, not the wavelet phasor mean field, and can be used to help attribute synchrony to particular causes. Thus the two techniques are often best used together: significance of phase synchrony can be identified with the wavelet phasor mean field, and then once significance is identified, synchrony can be described and studied using the wavelet mean field.

The power and plotphase methods also work on wmf objects.

Examples of the wavelet mean field applied to real data are in Sheppard et al. (2016) and Sheppard et al. (2019).

Fourier surrogates

This procedure can be done using surrog with surrtype="fft". Because only the phases of the Fourier transform are randomized, not the magnitudes, autocorrelation properties of the surrogate time series are the same as those of x(t).

This procedure can be done using surrog with surrtype=fft and with syncpres=TRUE (for synchrony-preserving surrogates) or with syncpres=FALSE (for independent surrogates). Autocorrelation properties of individual time series are preserved, as for the N = 1 case covered above. If synchrony-preserving surrogates are used, all cross-correlation properties between time series are also preserved, because cross spectra are unchanged by the joint phase randomization. Therefore synchrony is preserved.

Fourier surrogates tend to have normal marginal distributions (Schreiber and Schmitz 2000). Therefore, to ensure fair comparisons between statistical descriptors (such as coherences) of real and surrogate datasets, Fourier surrogates should only be applied to time series that themselves have approximately normal marginals. The Box-Cox transformations implemented in cleandat can help normalize data prior to analysis. If data are difficult to normalize, or as an alternative, the amplitude-adjusted Fourier transform surrogates method of the next section can be used instead.

Amplitude-adjusted Fourier surrogates

Amplitude-adjusted Fourier surrogates are described elsewhere (Schreiber and Schmitz 2000). Either synchrony preserving (syncpres=TRUE) of independent (syncpres=FALSE) amplitude-adjusted Fourier (AAFT) surrogates can be obtained from surrog using surrtype="aaft". AAFT surrogates can be applied to non-normal data, and return time series with exactly the same marginal distributions as the original time series. AAFT surrogates have approximately the same power spectral (and cross-spectral, in the case of syncpres=TRUE) properties as the original data.

Fast coherence

The fast coherence algorithm implemented in coh (option sigmethod="fast") implements Fourier surrogates only, and only applies for norm equal to none, powall, or powind. It is described in detail elsewhere (Sheppard, Reid, and Reuman 2017).

Alternatives to the “quick” method of assessing significance of wavelet phasor mean field values

When sigmethod is fft in a call to wpmf, the empirical wavelet phasor mean field is compared to wavelet phasor mean fields of Fourier surrogate datasets. The signif slot of the output is a list with first element "fft", second element equal to nrand, and third element the fraction of surrogate-based wavelet phasor mean field magnitudes that the empirical wavelet phasor mean field magnitude is greater than (a times by timescales matrix). For sigmethod equal to aaft, AAFT surrogates are used instead. Non-synchrony-preserving surrogates are used.

Model construction tools

First create a diver variable composed of an oscillation of period 12 years and an oscillation of period 3 years, and normally distributed white noise of mean 0 and standard deviation 1.5.

lts<-12
sts<-3
mats<-3
times<-seq(from=-mats,to=100)
ts1<-sin(2*pi*times/lts)
ts2<-sin(2*pi*times/sts)
numlocs<-10
d1<-matrix(NA,numlocs,length(times)) #the first driver
for (counter in 1:numlocs)
{
  d1[counter,]<-ts1+ts2+rnorm(length(times),mean=0,sd=1.5)
}

Next create a second driver, again composed of an oscillation of period 12 years and an oscillation of period 3 years, and normally distributed white noise of mean 0 and standard deviation 1.5.

ts1<-sin(2*pi*times/lts)
ts2<-sin(2*pi*times/sts)
d2<-matrix(NA,numlocs,length(times)) #the second driver
for (counter in 1:numlocs)
{
  d2[counter,]<-ts1+ts2+rnorm(length(times),mean=0,sd=1.5)
}

Next create an irrelevant environmental variable. With real data, of course, one will not necessarily know in advance whether an environmental variable is irrelevant to a population system. But, for the purpose of demonstrating the methods, we are playing the dual role of data creator and analyst.

dirrel<-matrix(NA,numlocs,length(times)) #the irrelevant env var
for (counter in 1:numlocs)
{
  dirrel[counter,]<-rnorm(length(times),mean=0,sd=1.5)
}

The population in each location is a combination of the two drivers, plus local variability. Driver 1 is averaged over 3 time steps in its influence on the populations, so only the period-12 variability in driver 1 influences the populations.

pops<-matrix(NA,numlocs,length(times)) #the populations
for (counter in (mats+1):length(times))
{
  aff1<-apply(FUN=mean,X=d1[,(counter-mats):(counter-1)],MARGIN=1)
  aff2<-d2[,counter-1]
  pops[,counter]<-aff1+aff2+rnorm(numlocs,mean=0,sd=3)
}
pops<-pops[,times>=0]
d1<-d1[,times>=0]
d2<-d2[,times>=0]
dirrel<-dirrel[,times>=0]
times<-times[times>=0]

If only the data were available and we were unaware of how they were generated, we may want to infer the causes of synchrony and its timescale-specific patterns in the populations. The wavelet mean fields of pops, d1 and d2 show some synchrony at timescales of about 3 and 12 for all three variables.

dat<-list(pops=pops,d1=d1,d2=d2,dirrel=dirrel)
dat<-lapply(FUN=function(x){cleandat(x,times,1)$cdat},X=dat)
wmfpop<-wmf(dat$pops,times,scale.max.input=28)
plotmag(wmfpop)

wmfd1<-wmf(dat$d1,times,scale.max.input=28)
plotmag(wmfd1)

wmfd2<-wmf(dat$d2,times,scale.max.input=28)
plotmag(wmfd2)

Thus we cannot know for sure from the wavelet mean fields whether population synchrony at each timescale is due to synchrony in d1, d2, or both drivers at that timescale. However, we can fit wavelet linear models.

Start by fitting a model with all three predictors. Only the "powall" option for norm is implemented so far.

wlm_all<-wlm(dat,times,resp=1,pred=2:4,norm="powall",scale.max.input=28)

We will carry out analyses for this model at long timescales (11 to 13 years) and short timescales (2 to 4 years) simultaneously. First test whether we can drop each variable.

wlm_all_dropi<-wlmtest(wlm_all,drop="dirrel",sigmethod="fft",nrand=100)
wlm_all_drop1<-wlmtest(wlm_all,drop="d1",sigmethod="fft",nrand=100)
wlm_all_drop2<-wlmtest(wlm_all,drop="d2",sigmethod="fft",nrand=100)

Examine results for dropping dirrel, long and short timescales. We find that dirrel does not need to be retained in either long- or short-timescale models, as expected given how data were constructed:

blong<-c(11,13)
bshort<-c(2,4)
wlm_all_dropi<-bandtest(wlm_all_dropi,band=blong)
wlm_all_dropi<-bandtest(wlm_all_dropi,band=bshort)
plotmag(wlm_all_dropi)

plotrank(wlm_all_dropi)

Examine results for dropping d1, long and short timescales. We find that d1 should be retained in a long-timescale model but need not be retained in a short-timescale model, again as expected:

wlm_all_drop1<-bandtest(wlm_all_drop1,band=blong)
wlm_all_drop1<-bandtest(wlm_all_drop1,band=bshort)
plotmag(wlm_all_drop1)

plotrank(wlm_all_drop1)

Examine results for dropping d2, long and short timescales. We find that d2 should be retained in both a short-timescale model and in a long-timescale model, again as expected:

wlm_all_drop2<-bandtest(wlm_all_drop2,band=blong)
wlm_all_drop2<-bandtest(wlm_all_drop2,band=bshort)
plotmag(wlm_all_drop2)

plotrank(wlm_all_drop2)

Note that only 100 randomizations were used in this example. This is for speed - in a real analysis, at least 1000 randomizations should typically be performed, and preferably at least 10000.

Amounts of synchrony explained

Now we have constructed models for short timescales (2 − 4 years) and long timescales (11 − 13 years) for the example, finding, as expected, that d1 is a driver at long timescales only and d2 is a driver at short and long timescales. How much of the synchrony in the response variable is explained by these drivers for each timescale band?

For short timescales, almost all the synchrony that can be explained is explained by Moran effects of d2:

se<-syncexpl(wlm_all)
se_short<-se[se$timescales>=bshort[1] & se$timescales<=bshort[2],]
round(100*colMeans(se_short[,c(3:12)])/mean(se_short$sync),4)

##     syncexpl   crossterms       resids           d1           d2       dirrel 
##      61.0564       7.0539      31.8896       0.9406      53.1484       0.2140 
## interactions        d1_d2    d1_dirrel    d2_dirrel 
##       6.7534       6.8652      -0.1210       0.0092

These are percentages of sychrony explained by various factors: syncexpl is total synchrony explained by the predictors for which we have data; crossterms must be small enough for the rest of the results to be interpretable; d1, d2 and dirrel are percentages of sychrony explained by those predictors; interactions is percentage of synchrony explained by interactions between predictors (see Sheppard et al. (2019)); and the remaining terms are percentages of synchrony explained by individual interactions.

For long timescales, Moran effects of both drivers are present, as are interactions between these Moran effects:

se_long<-se[se$timescales>=blong[1] & se$timescales<=blong[2],]
round(100*colMeans(se_long[,c(3:12)])/mean(se_long$sync),4)

##     syncexpl   crossterms       resids           d1           d2       dirrel 
##      94.4137       4.8170       0.7693      25.9920      29.5607       0.0104 
## interactions        d1_d2    d1_dirrel    d2_dirrel 
##      38.8506      38.5144       0.2230       0.1131

Note that cross terms are fairly small in both these analyses compared to synchrony explained. Results can only be interpreted when this is the case. See Sheppard et al. (2019) for detailed information on cross terms and interacting Moran effects.

The pattern of synchrony that would pertain if the only drivers of synchrony were those included in a model can also be produced, and compared to the actual pattern of synchrony (as represented by the wavelet mean field) to help evaluate the model.

pres<-predsync(wlm_all)
plotmag(pres)

plotmag(wmfpop)

The similarity is pretty good. Now make the comparison using the model with sole predictor d1.

wlm_d1<-wlm(dat,times,resp=1,pred=2,norm="powall",scale.max.input=28)
pres<-predsync(wlm_d1)
plotmag(pres)

The similarity with the wavelet mean field of the populations is pretty good at long timescales (where the model with sole predictor d1 was found to be a good model), but not at short timescales.

The synchrony matrix

There are numerous ways to generate a synchrony matrix, and synmat provides several alternatives. For an initial demonstration, create some data in two synchronous clusters.

N<-5
Tmax<-100
rho<-0.5
sig<-matrix(rho,N,N)
diag(sig)<-1
d<-t(cbind(mvtnorm::rmvnorm(Tmax,mean=rep(0,N),sigma=sig),
           mvtnorm::rmvnorm(Tmax,mean=rep(0,N),sigma=sig)))
d<-cleandat(d,1:Tmax,1)$cdat

Then make a synchrony matrix using Pearson correlation.

sm<-synmat(d,1:Tmax,method="pearson")
fields::image.plot(1:10,1:10,sm,col=heat.colors(20))

The function synmat provides many other options, beyond correlation, for different kinds of synchrony matrices. We demonstrate a frequency-specific approach. First create some artificial data.

N<-20
Tmax<-500
tim<-1:Tmax

ts1<-sin(2*pi*tim/5)
ts1s<-sin(2*pi*tim/5+pi/2)
ts2<-sin(2*pi*tim/12)
ts2s<-sin(2*pi*tim/12+pi/2)

gp1A<-1:5
gp1B<-6:10
gp2A<-11:15
gp2B<-16:20

d<-matrix(NA,Tmax,N)
d[,c(gp1A,gp1B)]<-ts1
d[,c(gp2A,gp2B)]<-ts1s
d[,c(gp1A,gp2A)]<-d[,c(gp1A,gp2A)]+matrix(ts2,Tmax,N/2)
d[,c(gp1B,gp2B)]<-d[,c(gp1B,gp2B)]+matrix(ts2s,Tmax,N/2)
d<-d+matrix(rnorm(Tmax*N,0,2),Tmax,N)
d<-t(d)

d<-cleandat(d,1:Tmax,1)$cdat

These data have period-5 oscillations which are synchronous within location groups 1 and 2, but are asynchronous between these groups. Superimposed on the period-5 oscillations are period-12 oscillations which are synchronous within location groups A and B, but are asynchronous between these groups. Groups 1 and 2 are locations 1 − 10 and 11 − 20, respectively. Group A is locations 1 − 5 and 11 − 15. Group B is locations 6 − 10 and 16 − 20. So the spatial structure of period-5 oscillations differs from that of period-12 oscillations. Strong local noise is superimposed on top of the periodic oscillations.

We measure synchrony matrices using portions of the the cross-wavelet transform centered on periods 5, and 12 (in separate synchrony matrices), to detect the different structures on different timescales.

sm5<-synmat(dat=d,times=1:Tmax,method="ReXWT",tsrange=c(4,6))
fields::image.plot(1:N,1:N,sm5,col=heat.colors(20))

sm12<-synmat(dat=d,times=1:Tmax,method="ReXWT",tsrange=c(11,13))
fields::image.plot(1:N,1:N,sm12,col=heat.colors(20))

This timescale-specific approach reveals the structure of the data better than a correlation approach.

sm<-synmat(dat=d,times=1:Tmax,method="pearson")
fields::image.plot(1:N,1:N,sm,col=heat.colors(20))

Several additional synchrony measures with which synmat can construct synchrony matrices are described in the documentation of the function.

Important note: synchrony matrices can have negative values for some of the methods provided by synmat. This is appropriate, since correlation and other measures of synchrony can be negative, but it complicates cluster detection (see next section).

Clustering

The function clust computes network modules/clusters and helps keep information about them organized. That function is also the generator function for the clust class. The class has print and summary and set and get methods (see the help file for clust_methods). The clustering algorithm used is a slight adaptation of that of Newman (2006) - see the next section for details. We illustrate the use of clust using the artificial data from the second example of the previous section.

#make some artificial coordinates for the geographic locations of where data were measured
coords<-data.frame(X=c(rep(1,10),rep(2,10)),Y=rep(c(1:5,7:11),times=2))

#create clusters based on the 5-year timescale range and map them using coords
cl5<-clust(dat=d,times=1:Tmax,coords=coords,method="ReXWT",tsrange=c(4,6))
get_clusters(cl5) #the first element of the list is always all 1s - prior to any splits

## [[1]]
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## [[2]]
##  [1] 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1

#call the mapper here
plotmap(cl5)

#plot mean time series for each module
plot(get_times(cl5)[1:100],get_mns(cl5)[[2]][1,1:100],type='l',col='red',
     ylim=range(get_mns(cl5)),xlab="Time step",ylab="Mean pop.")
lines(get_times(cl5)[1:100],get_mns(cl5)[[2]][2,1:100],type='l',col='green')
legend(x="topright",legend=c("mod1","mod2"),lty=c(1,1),col=c("red","green"))

#create wavelet mean fields for each module and plot
cl5<-addwmfs(cl5)
plotmag(get_wmfs(cl5)[[2]][[1]])

plotmag(get_wmfs(cl5)[[2]][[2]])

#create clusters based on the 12-year timescale range and map them using coords
cl12<-clust(dat=d,times=1:Tmax,coords=coords,method="ReXWT",tsrange=c(11,13))
cl12$clusters

## [[1]]
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## [[2]]
##  [1] 2 2 2 2 2 1 1 1 1 1 2 2 2 2 2 1 1 1 1 1

#call the mapper here
plotmap(cl12)

#plot mean time series for each module
plot(get_times(cl12)[1:100],get_mns(cl12)[[2]][1,1:100],type='l',col='red',
     ylim=range(get_mns(cl12)),xlab="Time step",ylab="Mean pop.")
lines(get_times(cl12)[1:100],get_mns(cl12)[[2]][2,1:100],type='l',col='green')
legend(x="topright",legend=c("mod1","mod2"),lty=c(1,1),col=c("red","green"))

#create wavelet mean fields for each module and plot
cl12<-addwmfs(cl12)
plotmag(get_wmfs(cl12)[[2]][[1]])

plotmag(get_wmfs(cl12)[[2]][[2]])

Color intensity on maps of clusters indicates the strength of contribution of a node to its module - details in the next section. The function addwpmfs is similar to the function addwmfs demonstrated above, but adds wavelet phasor mean field information to a clust object.

The clustering algorithm

The clustering algorithm used by clust is implemented in cluseigen, and is a generalization of the algorithm of Newman (2006). The algorithm makes use of the concept of modularity. Modularity was defined by Newman (2006) for any unweighted, undirected network paired with a partitioning of the nodes into modules. The modularity is then a single numeric score which is higher for better partitionings, i.e., for partitionings that more effectively group nodes that are more heavily connected and separate nodes that are less connected. The ideal goal would be to find the partitioning that maximizes the modularity score, but this is computationally infeasible for realistically large networks (Newman 2006). Instead, the algorithm of Newman (2006), which is computationally very efficient, provides a good partitioning that is not guaranteed to be optimal but is typically close to optimal (Newman 2006). The modularity itself can be computed rapidly, given a partitioning, with the function modularity in wsyn.

Modularity is defined in Newman (2006) for unweighted networks, and the definition generalizes straightforwardly to weighted networks for which all weights are non-negative. But synchrony matrices can represent weighted networks for which some weights are allowed to be negative. A generalization of modularity for this more general type of network was defined by Gomez, Jensen, and Arenas (2009), and the modularity function computes this generalization. Values are the same as the definition of Newman (2006) for non-negatively weighted networks. The algorithm implemented by cluseigen is a slight generalization of the original Newman (2006) algorithm that applies to weighted networks for which some edges are allowed to have negative weight.

We here describe the generalized algorithm, the validity of which was realized by Lei Zhao. Let w_ij be the adjacency matrix for a network, so the ijth entry of this matrix is the weight of the edge between nodes i and j, 0 if there is no edge. Let C_i be the community/module to which node i is assigned and let C_j be the same for node j. The original definition of modularity, for non-negative w_ij, is $$Q=\frac{1}{2w}\sum_{ij} \left(w_{ij}-\frac{w_i w_j}{2w}\right)\delta(C_i,C_j),$$ where w_i = ∑_jw_ij, $w=\frac{1}{2}\sum_i w_i$, and δ is the Kronecker delta function, equal to 1 when C_i = C_j and 0 otherwise.

For the case of partitioning into two clusters, Newman (2006) notes that $\delta(C_i,C_j)=\frac{1}{2}(s_i s_j +1)$, where we define s_i to be 1 if node i is in group 1 and −1 if it is in group 2. Then $$Q=\frac{1}{4w}\sum_{ij} \left(w_{ij}-\frac{w_i w_j}{2w}\right)(s_i s_j +1),$$ and it is easy to show this is $$Q=\frac{1}{4w}\sum_{ij} \left(w_{ij}-\frac{w_i w_j}{2w}\right)s_i s_j.$$ Defining a matrix B such that $B_{ij}=w_{ij}-\frac{w_i w_j}{2w}$, we have $$Q=\frac{1}{4w}\mathbf{s}^\tau \mathbf{B} \mathbf{s},$$ where τ is transpose and the bold quantities are the vectors/matrices composed of the indexed quantities denoted with the same symbol. Newman (2006) presents an argument that a good way to come close to optimizing Q over possible partitions into two modules is to choose s_i to be the same sign as the ith entry of the leading eigenvector of B, if the leading eigenvalue is postive (otherwise the algorithm halts with no splits, returning the trivial “decomposition” into one module). This gives the first split, according to the Newman algorithm, of the network into modules. Subsequent splits are handled in a similar but not identical way described by Newman (2006). (In particular it is incorrect to simply delete edges connecting the two modules from the first split and apply the algorithm to the resulting graphs.)

The generalized modularity of Gomez, Jensen, and Arenas (2009) is defined as follows. Let w_ij⁺ = max (0, w_ij) and w_ij⁻ = max (0, −w_ij) so that w_ij = w_ij⁺ − w_ij⁻. Let w_i⁺ = ∑_jw_ij⁺, w_i⁻ = ∑_jw_ij⁻, $w^+ = \frac{1}{2} \sum_i w_i^+$, and $w^- = \frac{1}{2} \sum_i w_i^-$. Then Gomez, Jensen, and Arenas (2009) justifies the definitions $$Q^+ = \frac{1}{2w^+} \sum_{ij}\left(w_{ij}^+-\frac{w_i^+ w_j^+}{2w^+}\right)\delta(C_i,C_j),$$ $$Q^- = \frac{1}{2w^-} \sum_{ij}\left(w_{ij}^- - \frac{w_i^- w_j^-}{2w^-}\right)\delta(C_i,C_j),$$ and $$Q=\frac{2w^+}{2w^+ + 2w^-}Q^+ - \frac{2w^-}{2w^+ + 2w^-}Q^-.$$ This is a generalization of the old definition of Q, in the sense that it reduces to that definition in the case of non-negatively weighted networks. Gomez, Jensen, and Arenas (2009) provides a probabilistic interpretation that generalizes the probabilisticc interpretation of Newman (2006).

It is straightforward to show $$Q=\frac{1}{2w^+ + 2w^-} \sum_{ij} \left(w_{ij}-\left(\frac{w_i^+ w_j^+}{2w^+} - \frac{w_i^- w_j^-}{2w^-}\right)\right)\delta(C_i,C_j).$$ Again considering an initial 2-module split and define s_i and s_j as previously, we have $$Q=\frac{1}{4w^+ + 4w^-} \sum_{ij} \left(w_{ij}-\left(\frac{w_i^+ w_j^+}{2w^+} - \frac{w_i^- w_j^-}{2w^-}\right)\right)s_i s_j,$$ which can be written in matrix form as $$Q=\frac{1}{4w^+ + 4w^-} \mathbf{s}^\tau \mathbf{E} \mathbf{s}.$$ Because the generalized modularity of Gomez, Jensen, and Arenas (2009) can be written in the same matrix format as the modularity expression of Newman (2006), the same eigenvector-based algorithm for finding a close-to-optimal value of the modularity can be used.

The quantity $$\frac{1}{2w^+ + 2w^-} \sum_{j} \left(w_{ij}-\left(\frac{w_i^+ w_j^+}{2w^+} - \frac{w_i^- w_j^-}{2w^-}\right)\right)\delta(C_i,C_j)$$ is the contribution of node i to the modularity. It is the extent to which node i is more connected to other nodes in its module than expected by chance (Newman 2006; Gomez, Jensen, and Arenas 2009), and can be interpreted as a strength of membership of a node in its module. The plotmap function (demonstrated above) has an option for coloring nodes according to this quantity.