Estimating and testing sequential causal effects based on alternative G-formula: an observational study of the influence of early diagnosis on survival of cardia cancer

Cancer diagnosis is part of a complex stochastic process, in which patients' personal and social characteristics influence the choice of diagnosing methods, diagnosing methods in turn influence the initial assessment of cancer stage, cancer stage in turn influences the choice of treating methods, and treating methods in turn influence cancer outcomes such as cancer survival. To evaluate the performance of diagnoses, one needs to estimate and test the sequential causal effect (SCE) under a specified regime of diagnoses and treatments in such a complex observational study, where the data-generating mechanism is unknown and modeling is needed for statistical inference. In this article, we introduce a method of statistical modeling to estimate and test SCEs under regimes of treatments (diagnoses and treatments in cancer diagnosis) in complex observational studies. By applying the alternative G-formula, we express the SCE in terms of the point effects of treatments in the sequence, so that the modeling can be conducted via the point effects in the framework of single-point causal inference. We illustrate our method by a medical example of cancer diagnosis with data from a Swedish prognosis study of cardia cancer.</p

Properties of Bu1 and Bu2 subgroups as separated by GMM clustering are similar to those of bursting subgroups as determined by the 500 Hz intraburst frequency criteria.

(A) Responses to 400 ms long synthetic stimuli in Bu1 units have shorter latency, higher peak firing rate, and more complete and rapid adaptation than responses in Bu2 units. (B) Bu1 unit maximum firing rate was sensitive to the rate of sound onset. (C) Bu1 unit VS was higher than Bu2 VS and peaked at intermediate SAM rates. (D) A majority of Bu1 units were synchronized at 16 Hz or higher SAM rate, in contrast with RS, FS, and Bu2 groups. (E) CI was highest for Bu1 units, indicating a tendency for spikes to occur at nearly the same time on each repetition of the vocalization stimuli. (F) This tendency is also reflected in the Bu1 group’s right shifted (toward temporal encoding) H versus q curve for decoding based on the Victor–Purpura spike distance metric. Data underlying this figure can be found in S2 Data. CI, correlation index; FS, fast-spiking; ISI, interspike interval; PSTH, peristimulus time histogram; RS, regular-spiking; SAM, sinusoidal amplitude modulation; VS, vector strength. (TIF)</p

Response properties to SAM, classified by GMM and all bursting units.

Same plots as Fig 6, but using unit type labels from the clustering analysis rather than the labels generated by criteria. Bursting units identified by 3 methods are shown for comparison: Bu (from GMM), PBu (from prestimulus logISIdrop), and Bucrit (from method of criteria, Bu1 and Bu2 combined). (A) Violin plot of maximum VS for RS (74), FS (65), Bu (40), PBu (45), and Bucrit (35) units. (B) Mean VS versus SAM modulation rate. (C) Violin plot of maximum synchronized rate for each unit type. (D) Fraction of responsive units that were synchronized at or above 4 Hz. (E) Fraction of responsive units that were synchronized at or above 16 Hz. (F) Average period histograms for stimulation at 2 Hz. Data underlying this figure can be found in S2 Data. Bu, bursting; FS, fast-spiking; GMM, Gaussian mixture model; RS, regular-spiking; SAM, sinusoidal amplitude modulation; VS, vector strength. (TIF)</p

Bursting units had the most strongly adapting responses to sustained stimuli at best frequency and were sensitive to the rate of sound onset.

(A) Smoothed PSTHs in response to 12.5, 50, and 400 ms stimuli at best frequency, bandwidth, and sound level. Bursting units had rapidly rising and strongly adapting onset responses. The effect is clearest for the Bu1 subgroup (intraburst frequency >500 Hz). Shaded bands shown standard errors; n values are RS (78), FS (44), Bu1 (12), and Bu2 (20). (B) An adaptation index was calculated from responses to 200 ms best frequency stimuli (larger values indicate more adaptation). Welch’s ANOVA revealed a significant group difference (F3,76.6 = 9.2, p 500 Hz) appear above the dashed line. (D) Driven rate, calculated over the entire response window and normalized to the maximum, peaks at 25 ms for Bu1 (due to lack of sustained response). In response to stimuli with increasingly slow onset ramps (E), Bu1 unit spiking became more distributed (example shown in (G)) or even nonresponsive. Maximum PSTH height decreased substantially with slower stimulus onsets for Bu1 units, but not for the other types (F). n values for ramp rate are RS (57), FS (29), Bu1 (11), and Bu2 (16). Recentered response heatmaps in (H), (I), (J), and (K) show temporal and receptive field differences between the unit types. Data underlying this figure can be found in S1 Data. PSTH, peristimulus time histogram; RS, regular-spiking.</p

Responses of bursting units to the “Mixed Vocalizations List” (S6A Fig).

Examples of diverse precise responses to vocalizations from bursting units (intraburst frequency shown in top right corner). Alternating light aqua shading indicates the stimulus duration. Data underlying this figure can be found in S2 Data. (TIF)</p

RS, FS, and Bu units had distinct types of responses to vocalizations analogous to their responses to sustained and SAM stimuli.

Responses are shown to “Call Type Lists” consisting of example tokens of the sustained phee call or the rapidly fluctuating trill call. Example spectrograms, waveforms, and corresponding response raster plots for bursting units are shown on the right, while complete stimulus spectrograms can be seen in S6B and S6C Fig. (A) Columns show responses of 4 example units of each type to 10 prerecorded phees. Bu units responded at particular moments, often at the onset, while RS and FS units were excited or inhibited throughout the call. (B) A Bu unit phase-locked to the rapid modulation (approximately 30 Hz) in some examples of trill calls. An FS unit was also capable of representing the modulation, but with many more spikes per cycle. Units shown were consistently labeled by criteria and GMM. The 5 distinct Bu units had intraburst frequencies of 476, 476, 526, 588, and 909 Hz (top to bottom). Alternating light aqua shading indicates stimulus duration. Bu, bursting; FS, fast-spiking; GMM, Gaussian mixture model; RS, regular-spiking; SAM, sinusoidal amplitude modulation.</p

In response to sinusoidally amplitude modulated sounds, bursting neurons had higher VS, were more likely to be synchronized and up to higher rates, and responded in the early phase of the cycle.

(A) Violin plot of maximum VS. Nonsignificant VS values were set to 0 (see Materials and methods). Bursting types had high VS values. For (A) and (C), Welch’s ANOVA was statistically significant (F3,42.0 = 11.3, p S1 Data. Bu, bursting; FS, fast-spiking; GMM, Gaussian mixture model; RS, regular-spiking; SAM, sinusoidal amplitude modulation; VS, vector strength.</p

Unit types did not differ grossly in terms of frequency, depth, or regional distribution despite showing consistent functional differences.

(A) Schematic showing location of auditory cortex along the lateral sulcus of the left hemisphere of the marmoset brain and cortical areas within the core region (dark shading) and belt region (light shading), based on (1–3). Within the core, the tonotopy experiences a low frequency reversal at the lateral border between AI and R, and a mid-frequency reversal at the medial border between R and RT. Recordings were made along the length of the lateral sulcus, primarily in core areas AI, R, and RT, with some likely inclusion of anterior and caudal belt. (B and C) BF and recording depth distributions were generally overlapping for the various unit types (RS, red circles; FS, blue squares; Bu1, dark green triangles; and Bu2 light green diamonds), and should not be a confounding cause for consistent unit type differences observed between unit types. Depths are expressed relative to the first spiking unit encountered from a superficial approach and were biased toward superficial layers due to the long recording times spent with each unit. One-way ANOVAs did not show a statistically significant difference in BF or depth between at least 2 groups (F(3,329) = 1.97, p = 0.12 and F(3,355) = 1.6, p = 0.19). (D andE) Maps of best frequencies of recorded units in the 2 marmosets used in this study, spanning from the low frequency region of anterior RT to the high frequency region of posterior AI. See (H and I) for scale. A small jitter was added to offset multiple units within the same track for visibility. Light gray x’s indicate units that could not be well driven by sound. (F and G) Unit types were distributed throughout recorded areas. For instance, Bu1 units (dark green triangles) were interleaved with other unit types. (H and I) When units were projected onto the sulcal axis, bursting units, and in particular Bu1 units, had higher CImax values regardless of anterior-posterior location. Data underlying this figure can be found in S2 Data. AI, primary auditory cortex; AL, anterolateral belt; BF, best frequency; CL, caudolateral belt; CM, caudomedial belt; FS, fast-spiking; ML, middle lateral belt; MM, middle medial belt; R, rostral core; RM, rostromedial belt; RS, regular-spiking; RT, rostrotemporal core; RTL, rostrotemporal-lateral belt; RTM, rostrotemporal-medial belt; (TIF)</p

S2 Data -

“S2_Data.xlsx” includes the data underlying the Supporting information figures. (XLSX)</p

Unsupervised classification also detects 3 major classes with strong agreement with criteria-based classification.

PCA was performed on a set of 8 features. The loading plot (A) revealed 2 main subsets of features, those pertaining to spike timing and burstiness, and those pertaining to fast spiking. We fit a 3-component GMM to the first 3 PCA components (B). The color of the points indicates the labels assigned by the previous method using criteria (triangular and diamond markers indicate Bu1 and Bu2 subgroups, respectively). The elongated ovoids represent the 3D contours of each of the 3 components of the GMM at half height. (D) The plot of the AIC and BIC versus number of components both support the choice of 3 components. (E) The data set was split into training and test sets (repeated 20 times), and the negative log-likelihood was calculated as a proxy for quality of fit. The negative log-likelihood increased after 3 components for the test set, suggesting that more than 3 components resulted in overfitting. (F) The confusion matrix between the type labels assigned by criteria and GMM showed high agreement. Data underlying this figure can be found in S1 Data. AIC, Akaike information criterion; BIC, Bayesian information criterion; GMM, Gaussian mixture model; PCA, principal component analysis.</p