Broad Spectrum epidemiological contribution of cannabis and other substances to the teratological profile of northern New South Wales: geospatial and causal inference analysis

Abstract

Background

Whilst cannabis commercialization is occurring rapidly guided by highly individualistic public narratives, evidence that all congenital anomalies (CA) increase alongside cannabis use in Canada, a link with 21 CA’s in Hawaii, and rising CA’s in Colorado indicate that transgenerational effects can be significant and impact public health. It was therefore important to study Northern New South Wales (NNSW) where cannabis use is high.

Methods

Design: Cohort. 2008–2015. Setting: NNSW and Queensland (QLD), Australia. Participants. Whole populations. Exposures. Tobacco, alcohol, cannabis. Source: National Drug Strategy Household Surveys 2010, 2013. Main Outcomes. CA Rates. NNSW-QLD comparisons. Geospatial and causal regression.

Results

Cardiovascular, respiratory and gastrointestinal anomalies rose with falling tobacco and alcohol but rising cannabis use rates across Queensland. Maternal age NNSW-QLD was not different (2008–2015: 4265/22084 v. 96,473/490514 > 35 years/total, Chi.Sq. = 1.687, P = 0.194). A higher rate of NNSW cannabis-related than cannabis-unrelated defects occurred (prevalence ratio (PR) = 2.13, 95%C.I. 1.80–2.52, P = 3.24 × 10^− 19). CA’s rose more potently with rising cannabis than with rising tobacco or alcohol use. Exomphalos and gastroschisis had the highest NNSW:QLD PR (6.29(2.94–13.48) and 5.85(3.54–9.67)) and attributable fraction in the exposed (84.11%(65.95–92.58%) and 82.91%(71.75–89.66%), P = 2.83 × 10^− 8 and P = 5.62 × 10^− 15). In multivariable geospatial models cannabis was significantly linked with cardiovascular (atrial septal defect, ventricular septal defect, tetralogy of Fallot, patent ductus arteriosus), genetic (chromosomal defects, Downs syndrome), gastrointestinal (small intestinal atresia), body wall (gastroschisis, diaphragmatic hernia) and other (hypospadias) (AVTPCDSGDH) CA’s. In linear modelling cannabis use was significantly linked with anal stenosis, congenital hydrocephalus and Turner syndrome (ACT) and was significantly linked in borderline significant models (model P < 0.1) with microtia, microphthalmia, and transposition of the great vessels. At robust and mixed effects inverse probability weighted multivariable regression cannabis was related to 18 defects. 16/17 E-Values in spatial models were > 1.25 ranging up to 5.2 × 10¹³ making uncontrolled confounding unlikely.

Conclusions

These results suggest that population level CA’s react more strongly to small rises in cannabis use than tobacco or alcohol; cardiovascular, chromosomal, body wall and gastrointestinal CA’s rise significantly with small increases in cannabis use; that cannabis is a bivariate correlate of AVTPCDSGDH and ACT anomalies, is robust to adjustment for other substances; and is causal.

Background

With major tobacco companies entering the cannabis market it is clear that cannabis commercialization is well under way [1]. Whilst much of the discussion relating to cannabis use and cannabis control is notably self-referential recent epidemiological reports suggest that intergenerational effects may be both significant and powerful enough to impact population-level health outcomes.

A recent report on Canada demonstrated that total congenital defects were three times more common in the northern Territories which smoked more cannabis than the Provinces and that the association was robust to socioeconomic adjustment [2]. A recent study from Colorado across the period of cannabis legalization showed that many defects rose parallel to increased cannabis consumption including all chromosomal defects (ACD), Downs syndrome and several cardiovascular defects including atrial septal defect (ASD) and patent ductus arteriosus (PDA), common defects which had not been previously linked with prenatal cannabis exposure (PCE) [3]. It was calculated that in Colorado over 11,000 extra defects occurred 2000–2014 related to increased cannabis use [3]. An Hawaiian study found that 21 defects were increased in mothers who were exposed only to cannabis [4].

Whilst some of these studies have used sophisticated geospatial modelling techniques [2] all epidemiological research is fundamentally associational in nature. However similar findings elsewhere strengthens the evidence base.

Northern New South Wales (NNSW, NSW) is a well known drug using and cannabis cultivation area 760 km from Sydney but only 180 km from tertiary pediatric care centres in Brisbane and 111 km from Southport both in Queensland (QLD). Although lying within New South Wales administratively many of its neonatal CA’s are evacuated to tertiary pediatric hospitals in Queensland under the Neonatal Retrieval Scheme (NRS) [5] and their data thus appears in Queensland statistics. This therefore presents an ideal opportunity to directly compare NNSW and Queensland neonatal epidemiology.

Our hypothesis was that cannabis use would be associated with increased congenital anomalies and was formulated prior to study commencement.

Methods

Data

Data on congenital anomaly rates for Queensland Health service areas including northern New South Wales was taken from the Congenital Anomaly Linked File (CALF) from Queensland Health [6]. Annual data by area has not been publicly released. Data on maternal age was from the QLD and NSW annual Mothers and Babies reports [7, 8]. CALF data includes numbers, rates and confidence intervals for the data. Drug use data for last month cigarette use, last month binge alcohol and last year cannabis use by area was obtained from the Australian Institute of Health and Welfare from the National Drug Strategy Household Survey (NDSHS) 2010 and 2013 [9] and averaged to obtain a mean rate by area across this period pursuant to our custom data request. Data was matched manually between drug use and congenital anomaly datasets. Areal shapefiles were taken from the Australian Government national website [10]. The northern coastal area of NSW was added on to the Queensland Health shapefile. This depiction of the NNSW catchment area is illustrative only and not intended to be exact as the geographic boundaries of the NRS are not defined [5].

Congenital anomalies were defined as cannabis related or not based on a literature review and recent reports [2, 3, 11,12,13] particularly [4].

Sample characteristics

The samples were whole population samples and included all births in all health regions of Queensland. Hence inclusion and exclusion criteria were not applied as the samples were complete population samples. Sample capture in Northern New South Wales (NSW) appears to have been incomplete as some birth defects were not captured in the Queensland data and were managed in New South Wales. Hence the Northern New South Wales rates described in the present report clearly represents underestimates of the total rates. As the NSW and Queensland datasets are not directly comparable it is not possible to directly merge the two datasets to form a complete picture. This study limitation is discussed further in the limitations paragraph of the Discussion.

Statistics

Data was processed in RStudio version 1.2.1335 based on R version 3.6.1 on 16th April 2020. Two-by-two tables were analyzed in package epiR using epi.2by2. Graphs were drawn with R-Base and ggplot2 [14] and in Microsoft Excel. Maps were drawn using the R packages ggplot2 and sf (“simple features”) [14, 15]. The software is freely available online and is directly loaded from inside RStudio software downloaded from the Comprehensive “R” Archive Network (CRAN). Licences for all R software are freely available with the R packages sourced from CRAN and its various mirrors internationally. Principal component analysis was conducted using the psych package. Linear regression was performed in Base-R. Batch extraction of all linear model coefficients by different defects was performed with broom and purrr packages. Correction of P-values for multiple testing was applied using the corrections of Holm, Bonferroni, Benjamini and Yekutieli, False Discovery Rate and Hommel as indicated. Links between neighbouring areas sharing an edge or corner (“queen”-relationships by analogy with chess moves) were derived with the poly2nb function from spdep and edited as indicated. This neighbourhood map was used to calculate the geospatial weights matrix for spatial regression.

Geospatial regression was performed with the spreml function from package splm [13, 16] using the derived spatial weights matrix. Spatial models may include the parameters phi, psi, rho and lambda for random effects in the error term, serial autocorrelation in the residuals, spatial errors, and autocorrelation in the spatial errors respectively. All spatial models used a full error structure of Kapoor Kelejian and Prucha [17] and had serially correlated remainder errors and random effects (sem2srre). The appropriateness of this error structure was formally tested by substituting various alternative forms and comparing results including the maximum likelihood ratio at model optimization (LogLik in the tables) and spatial Hausman tests and by examining the significance of the final model parameters [18]. Models were spatially lagged and not lagged as indicated.

Inverse probability weights (IPW) were derived with the ipw package in R using cannabis use as the exposure of interest, tobacco in the numerator and tobacco and alcohol in the denominator. IPW weights were then used in robust regression models conducted in the R package survey, and in mixed effects models in the R package nlme to generate datasets pseudo-randomized for cannabis exposure. This allowed causal relationships to be assessed. E-Values were calculated with the R package EValue to quantify the degree of association some unmeasured confounder would require with both dependent and independent variables to explain away the observed effect.

It was necessary to use various regression model structures for several reasons. Only spatial models can capture the spatial effects amongst the data, however they cannot be inverse probability weighted, so that causality can not be studied directly from them. Spatial models do provide a model variance which allows e-Values to be calculated from them. Robust marginal structural models can be inverse probability weighted, but do not provide model standard deviations so that e-Values cannot be calculated from them directly. Mixed effects models can be inverse probability weighted and also provide model standard deviations, but do not capture spatial effects and are also not robust in design. Hence for a full picture and understanding of the causal and spatial processes involved one needs to consider output from all of these various model structures together with their applicable e-Values.

For all regression models model reduction was by the classical method with sequential deletion of the least significant term. Missing data was casewise deleted at multivariable regression. P < 0.05 was considered significant.

Data availability statement

Data including R software code have been made available online in the Mendeley Data Repository at https://doi.org/10.17632/cjzfyktz5m.1.

Ethics

This study was approved by the Human Research Ethics Committee of the University of Western Australia on 15th April 2020 (No. RA/4/20/4724).

Results

Input data is shown in an online supplementary csv file. Supplementary Table 1 provides comparative congenital anomaly data between QLD and NSW by both numbers and rates including defect relationship to cannabis [6]. Denominator data was calculated from the numbers and rates supplied in that file. It was verified from the annual QLD Health Mothers and Babies reports 2008–2015 which show 509,095 births in this period [8]. The “Interstate and Overseas” designation in the CALF file includes offshore islands such as Christmas, Norfolk, Cocos and Lord Howe Islands which together have a population of 4518. The prime catchment area of the NRS is Northern NSW which has a population of 296,531 [10]. Hence only 1.5% of the population in this designation is likely to come from outside NNSW. The view that the “Interstate and Overseas” designation refers primarily to NNSW is confirmed by QLD Health Ministerial correspondence (Minister Steven Miles, 05/04/2018). The denominator figure calculated for NNSW in this manner is 4800 births.

It should be noted that NNSW birth defect data also appears in NSW Health records [7]. One notes that the rates of congenital anomalies reported for this region in the NSW Mothers and Babies reports are about half those of the rest of the state. This is presumably related to the relocation of many cases into Queensland through the NRS. Queensland congenital anomaly rates are much higher than those reported elsewhere so it is not possible simply to combine NSW and QLD Health reports. Therefore this report is limited to consideration of the QLD Health CALF file only.

Drug use data is shown in Supplementary Table 2. It is noteworthy that the Richmond-Tweed NNSW area has a middle ranking for tobacco and alcohol use, but a first ranking for cannabis use.

Maternal age is a major factor bearing on congenital anomaly rates and it is known to be strongly linked with chromosomal anomaly rates. For years 2008–2015 4265/22084 (19.31%) births in NNSW were to mothers > 35 years compared to 96,493/490514 (19.67%; 18,581 missing maternal age data) in Queensland (Chi.Sq. = 1.687, df = 1, P = 0.194). This compares to 22,133/92242 (23.9%) of 2012 births in the rest of NSW (Chi.Sq. = 98.954, df = 1, P = 4.33 × 10^− 10) indicating that NNSW mothers are younger in than elsewhere in NSW.

Interestingly CALF Table 1 shows rises in the rates of several defects (Fig. 1) including CVS defects, atrial septal defect (ASD) and ventricular septal defect (VSD) which are highly significant (Supplementary Table 3). Intriguingly the mean incidence of daily smoking tobacco and high risk alcohol use dropped across this period and annual cannabis use rose from 10.5 to 11.3%. The principal component of the combination of cardiovascular, gastrointestinal and respiratory anomalies also rose significantly across this period. These data suggest that cannabis may be a more potent and more important teratogen than tobacco and alcohol.

Table 1 Results of geospatial additive regression by selected congenital anomalies - spatial error models

Fig. 1

Trend analysis over time of (a) drugs and (b) various congenital defects and defect classes

Figure 2 shows a qualitative choropleth map-graph for the major CA classes. The yellow zones reflect high incidence and dark blue low incidence.

Fig. 2

Choropleth maps of congenital anomaly class rates across QLD and NNSW. High rates are shown in yellow and low rates in dark blue. Maps were drawn using R package “sf” [15]

Supplementary Figs. 1–3 present choropleth maps of CA incidence by area. Figure 3 shows chromosomal anomaly incidence for which data is available.

Fig. 3

Choropleth maps of various chromosomal anomaly rates across QLD and NNSW. High rates are shown in yellow and low rates in dark blue. Maps were drawn using R package “sf” [15]

Figure 4 was drawn in Excel and shows the confidence intervals from CALF for common, intermediate frequency and rare defects for cannabis-related (CRD) and cannabis not related (CNRD) defects. For most of the cannabis-unrelated defects the confidence intervals overlap. For most of the cannabis-related defects the confidence intervals either do not overlap, or are near the lower end of the QLD C.I.’s. Figure 5 expands this list for rare congenital anomalies and in general terms tends to continue this trend of non-overlapping confidence intervals or confidence intervals at the more extreme end of the range of the confidence intervals for cannabis-related defects.

Fig. 4

Confidence Intervals of cannabis-related and cannabis-unrelated congenital anomalies for (a) common, (b) intermediate frequency and (c) rare congenital anomalies

Fig. 5

Confidence Intervals of rare cannabis-related and cannabis-unrelated congenital anomalies in more detail

Figure 6 compares the QLD and NNSW CA rates based on the relative rates quoted in the Queensland CALF file.

Fig. 6

Relative Rates of Congenital Anomalies in Queensland v NNSW

Figure 7 compares all the rate ratios of defects using the quoted rates in the CALF file. Figure 8 makes a similar comparison with log rates and shows clearly that most of the cannabis-related defects are more common in NNSW.

Fig. 7

Relative Rate Ratios between Queensland and NNSW by relationship to Cannabis

Fig. 8

Log (Relative Rate Ratios) between Queensland and NNSW by relationship to Cannabis

CNRD were more common in QLD (23,737/509095, 4.66% v 185/4800, 3.85%, Chi Sq. = 7.002, df = 1, P = 0.0081). CRD were more common in NNSW (394/4800 8.21% v 16,346/509095, Chi Squ. = 376.86, df = 1, P = 6.01 × 10^− 84). CRD were more common in NNSW than CNRD (394/4800 v 185/4800, Prevalence Ratio (PR) = 2.13 (95%C.I. 1.80–2.52), Chi Sq. = 80.284, P = 3.24 × 10^− 19).

Supplementary Table 4 lists the PR’s, attributable fraction in the exposed (AFE) and attributable fraction in the population (AFP) along with their C.I.’s and applicable P-values for all defects and defect classes. They decline from exomphalos (PR = 6.29 (2.94–13.48), AFE = 84.11% (65.95–92.58%) and AFP = 4.71% (0.55–8.69%)) and gastroschisis (PR = 5.85 (3.54–9.67), AFE = 82.91% (71.75–89.66%) and AFP = 4.34% (1.79–6.82%)).

Figure 9 illustrates the PR’s and C.I.’s for CRD and CNRD. Figure 10 shows the AFE’s and C.I.’s for CRD and CNRD. Figure 11 shows the AFP’s and C.I.’s for CRD and CNRD. P-values are illustrated in Fig. 12.

Fig. 9

Prevalence Ratios (a) cannabis-related and (b) cannabis-unrelated congenital anomalies

Fig. 10

Attributable Fraction in the Exposed for (a) cannabis-related and (b) cannabis-unrelated congenital anomalies

Fig. 11

Attributable Fraction in the Population for (a) cannabis-related and (b) cannabis-unrelated congenital anomalies

Fig. 12

P-Values for (a) cannabis-related and (b) cannabis-unrelated congenital anomalies

Figure 13 shows five main defect classes charted against the use of tobacco, binge alcohol and cannabis. Rising trends with cannabis seem to apply to CNS, cardiovascular and chromosomal anomalies.

Fig. 13

Response of five major congenital anomaly classes to (a) Tobacco, (b) Alcohol and (c) cannabis exposure

Figure 14 charts all 55 anomalies and anomaly classes against tobacco use. Figure 15 performs a similar function for binge alcohol.

Fig. 14

Congenital Anomaly Rate by Cannabis Use Rate – Ordered by Slopes of Least Squares Regression Lines

Fig. 15

Congenital Anomaly Rate by Tobacco Use Rate – Ordered by Slopes of Least Squares Regression Lines

When a similar exercise is undertaken for cannabis exposure rising trends appear in several defects in the top three rows especially in cardiovasculature, chromosomal anomalies and body wall defects (Fig. 16).

Fig. 16

Congenital Anomaly Rate by Alcohol Use Rate – Ordered by Slopes of Least Squares Regression Lines

Supplementary Table 5 lists the regression coefficients and their significance levels in ascending order of P-values for cannabis exposure. Supplementary Tables 6 and 7 show this table listed in order of ascending P-values for tobacco and alcohol respectively.

Supplementary Table 8 lists the significant output of a linear regression where the defect rate was related to additive terms of tobacco, binge alcohol use and cannabis use and includes all P < 0.3. This procedure selects 18 defects for further study. For biological and epidemiological reasons Trisomies 13 and 18 were also included.

Spatial analysis algorithms do not tolerate missing data. Hence linear regression was used to investigate 8 defects where spatial data was incomplete. Supplementary Table 8 shows the results of a model interactive in substances. Cannabis use was identified as being linked with several defects in this table including Turners syndrome.

Figure 17 shows the geospatial relationships which were derived from spdep::poly2nb and then edited to include all geospatial links.

Fig. 17

Geospatial Links Between regions – (a) Edited and (b) Final. Maps were drawn using R package “sf” [15]

Table 1 gives the results of geospatial regression for a model with additive terms in drug exposure. These are spatial error models and are not spatially lagged. In the additive model series cannabis is independently linked with all eight anomalies particularly cardiovascular (ASD, PDA and tetralogy of Fallot, ToF) and chromosomal (ACD and Downs syndrome), gastroschisis and small intestinal atresia.

Table 2 shows the results of an interactive spatial model and finds that cannabis is more strongly linked with these same defects. VSD is now positively associated as is diaphragmatic hernia which have both been previously noted to be cannabis-associated [19, 20].

Table 2 Results of geospatial interactive regression by selected congenital anomalies spatial error models

A similar exercise is executed for spatially lagged spatial error (spatially autocorrelated with autocorrelated error components, SARAR) additive (Table 3) and interactive (Table 4) models with very similar results. In each case spatial error models were superior to combined SARAR models, as judged by the log maximum likelihood values, spatial Hausman tests and the largely non-significant lambda coefficients.

Table 3 Results of geospatial additive regression by selected congenital anomalies combined spatial error and spatial lag (SARAR) models

Table 4 Results of geospatial interactive regression by selected congenital anomalies combined spatial error and spatial lag (SARAR) models

One notes also that in a number of spatial error models spatial model parameters rho and lambda are noted to be highly significant. This therefore justifies the use of spatial models and also suggests that spatial factors are significant in considering clinical teratological patterns.

Table 5 summarizes the results of the above spatial models to facilitate comparisons between the various substances. The four spatial model structures are listed across the top of the Table and the various substances are listed in the rows. One notes that cannabis was independently predictive in 27 of the 36 models compared to only five each for tobacco and alcohol. Cannabis was involved in 19 interactive terms compared to ten and 13 for alcohol and tobacco respectively. It follows therefore that cannabis was involved in 46 terms compared to 15 and 18 for alcohol and tobacco respectively.

Table 5 Summary of geospatial regression model results by substance

Having demonstrated a strong associational relationship between drug exposure and several congenital anomalies the next issue of importance relates to the issue of whether the relationship was causal or not. Inverse probability weights were generated and used to derive a dataset pseudorandomized for cannabis exposure. Data was processed by robust interactive generalized linear modelling functions. As shown in Table 6 cannabis was significantly related to 18 anomalies either alone or in interaction with tobacco and alcohol.

Table 6 Robust generalized linear regression models

This exercise was repeated with (non-robust) mixed effects modelling as such models in R have standard deviations associated with them, which is required in the E-values algorithm which follows. As shown in Table 7 similar results were obtained for 11 congenital anomalies.

Table 7 Mixed effects model results for selected congenital anomalies

It is conceivable that the described relationships were related to some factor other than the measured covariates. E-Values quantitate the degree of association required of some unmeasured confounder with both cannabis exposure and the dependent variables to explain away the described effect. One notes that 1.25 is the value quoted in the literature as the threshold of interest in sensitivity analyses [21]. Table 8 lists 27 e-Values and finds that 21 minimum e-Values were greater than 1.25 and ranged up to 3.8 × 10³⁰ for geospatial models and up to infinity for mixed effects models, making uncontrolled confounding unlikely.

Table 8 E-values for key regression parameters

Discussion

This investigation presents many intriguing findings. Despite the several technical shortcomings of this dataset it is fascinating for the details and tantalizing clues which have been revealed. Importantly most of its major findings have been confirmed previously in other locations particularly in Colorado, Hawaii, Canada, and USA and by professional bodies such as AHA, AAP and CDC lending support to the strength of its principal results [2, 3, 11, 13, 19, 20].

NNSW has higher prevalence rates of the cannabis related anomalies: neural tube defects; small intestinal atresia; body wall defects: exomphalos, gastroschisis, diaphragmatic hernia; the cardiovascular disorders: ASD, VSD, PDA, tetralogy of Fallot, and transposition of the great vessels (TxGrVess); and the genetic disorders: all chromosomal disorders, Downs syndrome, Turners syndrome and trisomy 18. Amongst the defect classes cardiovascular, respiratory, and chromosomal anomalies were elevated. Some of these associations have been previously reported [3, 4, 22] and were seen in our unpublished analyses of US data.

QLD Health data showed that the NNSW CI’s for CRD’s were mostly non-overlapping or were at the extreme end of the QLD CI’s. CRD’s had higher rate ratios than CNRD’s.

Rising rates of cardiovascular, gastrointestinal and respiratory defects, and their first principal component were associated with falling rates of tobacco and alcohol use but rising cannabis use, just as was found in Colorado and USA [3].

At geospatial and linear regression the cardiovascular defects ASD, VSD, PDA, ToF, TxGrVess; the chromosomal defects ACD, Downs, Turners, Trisomy 13; the body wall defects gastroschisis, exomphalos, diaphragmatic hernia; the GI disorders small intestinal atresia and anal stenosis were all linked with cannabis exposure and for most cannabis exposure was an independent risk factor.

Rising rates of cannabis exposure were more strongly associated with cardiovascular, chromosomal, gastrointestinal and body wall defects than were rising rates of tobacco or alcohol exposure.

Analysis of this dataset by the formal techniques of causal inference analysis including inverse probability weighting and E-Values demonstrated that the described relationships fulfil the criteria for causal relationships.

These results show a striking concordance with epidemiological series from elsewhere. ASD, VSD, ToF, obstructive urinary disorders, hydrocephalus, anal anomalies and Downs syndrome were linked with PCE in a large Hawaiian series [4]. VSD has previously been linked with PCE [19]. Neural tube defects were noted to be elevated in a cannabis-related manner in Canada and Hawaii [4, 11]. ASD, PDA ACD and Downs were seen to rise in close temporal association with increased cannabis use in Colorado [3]. Exomphalos was implicated in animals [23, 24] and in some clinical series including in Queensland [25]. Transposition of the great vessels has previously been linked with paternal PCE [26]. Indeed in Canada total CA’s were linked with increased cannabis use after controlling for income and sociodemographic variables [2].

Many series implicate PCE in gastroschisis aetiology with a meta-analyzed bivariate O.R. = 4.12 (95%C.I. 3.45–4.91) [4, 27,28,29,30,31,32]. Our findings PR = 5.85 (3.54–9.67) contradict those of a 2011 NSW Health report on gastroschisis in this region [33] which erroneously applied an inflated Bonferroni correction to obviate a significant result. Indeed if the 9 cases reported in NSW [33] are added to the 16 cases reported in QLD the PR rises further to 9.13 (6.07–13.72).

Increasing reports from diverse sources indicate that the evidence is building that cannabis has significant teratological activities in humans in agreement with animal studies where many severe defects including oedema, exomphalos, phocomelia, spina bifida, myelocoele, exencephaly and foetal loss were documented [23, 24]. Concordant reports from Hawaii, Colorado and Canada suggest that the findings reported herein are indeed valid and are generalizable elsewhere. Given that likely half the NNSW congenital anomalies are reported internally within NSW [7] this suggests that the teratological situation in NNSW is indeed serious. Moreover some of the CA described here, especially chromosomal defects, are heavily therapeutically aborted antenatally again suggesting that the situation may well be much worse than our description suggests. Our analysis strongly implicates cannabis use as a likely underlying factor.

When one also considers the known epigenetic actions of cannabis [2, 12, 34,35,36,37] and its associations with developmental neurological dysfunction and autism [38,39,40,41,42] concerns relating to the intergenerational actions of cannabis are heightened.

From both the present data and from similar international analyses a number of important clinical implications arise. Notwithstanding its popular relatively benign image such analyses indicated not only that the potential teratological impacts of cannabis are significant but that they are likely causal in nature. Patients considering commencing a family should be encouraged to desist from all drugs prior to conception including cannabis. Patients who do fall pregnant and who are consuming cannabis should be encouraged to reduce and cease. Patients wishing to access treatment to assist with such withdrawal should be provided every encouragement and assistance to do so. Patients should be warned that the evidence base for the use of cannabis for most of its touted clinical indications is weak. Patients should be advised to avoid cannabis for morning sickness of pregnancy. Heavy cannabis smokers should be warned that cannabis hyperemesis can mimic hyperemesis gravidarum.

Moreover since the debate relating to cannabis is typically highly individualistic it seems prudent that medical professional organizations should partner with public health agencies and community groups to enlarge the focus of popular debate from the simply self-referential to a broader multigenerational perspective.

One major toxicological conclusion which follows directly from these studies is that access to cannabis should be highly restricted. Indeed such work calls into question the whole issue of the long term advisability of cannabis medicalization / legalization and the sustainability of such paradigms from a teratological perspective.

The present work has not considered neurological sequalae in the newborn and childhood as has previously been reported to overlap the autistic spectrum disorder and ADHD and thereby potentially play a major role in the modern widespread epidemic of these disorders [38,39,40, 43]. When such data is factored into consideration the imperatives for reconsideration and re-evaluation of cannabis legalization overall are largely increased.

Comparison with alcohol

It is of interest to summarize and compare some of these results directly between cannabis and alcohol as the latter is a known human foetal teratogen and many learned bodies recommend strongly against tobacco exposure in pregnancy [44].

Figure 13 is a scatterplot of the frequency of the defect classes compared to the various substances. This Figure shows clearly that increasing cannabis use is associated much more strongly with several classes of congenital anomalies in this dataset than either tobacco or alcohol. The regression lines in this figure slope upwards much more strongly for cannabis for CNS defects, cardiovascular defects and chromosomal defects than for either of the other two substances.

Figures 14, 15 and 16 perform a similar role for each individual defect by the three substances tobacco, alcohol and cannabis. On the tobacco and alcohol scatterplots most of the regression lines are flat or falling. In contradistinction on the cannabis scatterplot many of the first 32 defects appear to be rising with positive slope. This is quantified in Supplementary Tables 5–7 for cannabis, tobacco and alcohol respectively. The slopes of the first 8 CA’s are significant and seven slopes are positive for cannabis. This compared to tobacco and alcohol where only two and four slopes are significant respectively and all the slopes are negative.

Table 1 presents the remarkable result that of eight additive spatial models cannabis is independently predictive for all eight defects and indeed tobacco and alcohol do not appear in final models. Similarly in Table 3 cannabis is independently predictive for eight of nine defects in additive SARAR spatial models. Alcohol only features in the model for Downs syndrome and its regression coefficient is negative. These differences are compared directly in Table 5, where as noted cannabis is implicated in 46 terms compared to 15 for alcohol and 18 for tobacco. Cannabis is implicated independently in 27 terms compared to five each for tobacco and alcohol.

The overall conclusion then from this detailed comparison must be that cannabis is a relatively more powerful or more potent human teratogen than alcohol.

Causal inference

A classical criticism of correlative studies is that “correlation does not equal causation.” Judea Pearl, one of the leading causal statisticians in the world, has described this criticism as arising from what has been historically the “causalophobic” science of statistics [45, 46]. In relation to the present study the following points should be mentioned. Firstly to observe that an exposure and an outcome are associated not only statistically but also across space carries more weight than a simple statistical association. Secondly inverse probability weighting has been used in mixed effects and robust structural marginal models with very highly significant results. Inverse probability weighting is well established in the literature as transforming an observational study into a pseudo-randomized population from which casual inferences can properly be drawn. Thirdly we have used e-Values to quantify unmeasured confounding as a notorious source of extraneous confounding not controlled by the small number of covariates employed in the present analysis. E-Values provide a quantitative estimate of the degree of association required of any extraneous factor with both the exposure and the outcome to explain away the observed effect. Whilst in the literature e-Values above 1.25 have been stated to be noteworthy [21] our minimum e-Values ranged up to infinity in mixed effects models, and up to 5.2 × 10^− 13 in spatial models. This finding implies both the causal nature of the relationship, and also that the inclusion of further parameters in the model would not obviate the described effects.

Hence our study demonstrates a causal relationship of drug and particularly cannabis exposure to several congenital anomalies. The causal relationship in this case is greatly strengthened by the existence of similar results from other places in the world as described [2,3,4, 11, 19] and the existence of a plethora of biological and epigenetic processes to account for these effects as mentioned [3, 11,12,13, 34, 35, 37, 39, 47,48,49,50].

It is further noted that the present findings fulfill all of the qualitative and quantitative Hill criteria for causality [51].

Our study has several strengths and limitations. Its strengths include access to whole population data for Queensland and a significant portion of the NNSW data. The CA rates and confidence intervals were already provided by QLD Health. The NDSHS is a nationally representative survey conducted every three years and the authoritative source for most Australian drug use data. Our analytical strategy combined CA with drug exposure data which is unusual and useful. We have employed a variety of powerful statistical techniques in this investigation including geospatial analysis, inverse probability weighting, mixed models and E-Values. Study limitations relate mainly to the remote location of the NNSW area close to the Queensland border and the small numbers of some anomalies reported. Losses due to treatment within NSW and to stillbirths and prenatal therapeutic abortion occurring preferentially in CA babies implies that the present findings are conservative estimates. The very high CA rate reported in Queensland has not been explained despite formal enquiry. The origin of the NNSW denominator figure is unclear. NSW Mothers and Babies reports [7] indicate that during 2008–2015 22,084 babies were born in Northern NSW and 30,848 in the central coast region, totalling 52,932 births. These regions are shown together in our maps. Hence over 11 times the data is available as was used in this analysis if it can be properly collated between the two jurisdictions of NSW and Queensland. This would then facilitate geotemporospatial statistical modelling. This proper collation and assembly of data is a top research priority for future studies. The remote location of NNSW together with its somewhat trans-jurisdictional status has apparently made such a collation difficult in the past.

Conclusions

In conclusion study data indicate that prenatal cannabis exposure is a significant and robust covariate of many congenital anomalies in NNSW particularly affecting the cardiovascular, chromosomal, body wall and gastrointestinal systems and is highly significant for 10 cannabis-related defects on geospatial analysis. Close concordance between these results and previous reports from Hawaii, Colorado, and Canada and with unpublished USA studies suggest our findings are reliable and generalizable. On all studies cannabis teratogenesis seems to be more concerning than the established teratogens tobacco and alcohol. Fulfillment of the criteria for causal relationships has been demonstrated. Further geospatial epidemiological and basic science research is a priority given cannabis commercialization. Even beyond the obvious jurisdictional health cost-shifting implications careful and thorough further investigation of the teratological profile of NNSW by coordinated investigations between NSW and Queensland over time to current would appear to be a major international research priority with implications far beyond our shores.