Effects of cannabis on congenital limb anomalies in 14 European nations: A geospatiotemporal and causal inferential study
Division of Psychiatry, University of Western Australia, Crawley, WA 6009, Australia
School of Medical and Health Sciences, Edith Cowan University, Joondalup, WA 6027, Australia
*Correspondence address. 39 Gladstone Rd., Highgate Hill, Brisbane, Queensland, Australia. Tel: (+617) 3844-4000; Fax: (+617) 3844-4015; E-mail: stuart.reece@bigpond.comAbstract
Cannabinoid exposure is increasing in some European nations. Europe therefore provides an interesting test environment for the recently reported link between cannabis exposure and congenital limb anomaly (CLA) rates (CLARs). Exponential genotoxic dose–response relationships make this investigation both intriguing and imperative. Annual CLAR in 14 nations were from Epidemiological Surveillance of Congenital Anomalies. Drug use rates were from European Monitoring Centre for Drugs and Drug Dependency. Median household income was from the World Bank. E-values provide a quantitative measure of robustness of results to confounding by extraneous covariates. Inverse probability weighting is an important technique for equalizing exposures across countries and removing sources of bias. Rates of CLA, hip dysplasia and the whole group of limb anomalies were higher in countries with increasing daily cannabis use (P = 1.81 × 10−16, 0.0005 and 2.53 × 10−6, respectively). In additive inverse-probability-weighted panel models, the limb reduction-resin Δ9-tetrahydrocannabinol (THC) concentration E-value estimate was 519.93 [95% lower bound (mEV) 49.56], order Resin > Herb ≫ Tobacco > Alcohol. Elevations were noted in 86% E-value estimates and 70.2% of mEVs from 57 E-value pairs from inverse-probability-weighted panel models and from spatial models. As judged by the mEV the degree of association with metrics of cannabis exposure was hip dysplasia > polydactyly > syndactyly > limb anomalies > limb reductions with median E-value estimates from 3.40 × 1065 to 7.06 and median mEVs from 6.14 × 1033 to 3.41. Daily cannabis use interpolated was a more powerful metric of cannabis exposure than herb or resin THC exposure. Data indicate that metrics of cannabis exposure are closely linked with CLAR and satisfy epidemiological criteria for causality. Along with Hawaii and the USA, Europe now forms the third international population in which this causal link has been demonstrated. Cannabis as a predictor of limb anomalies was more potent than tobacco or alcohol. Cannabinoid access should be restricted to protect public health and the community genome/epigenome transgenerationally.
Introduction
Congenital anomalies (CAs) of limb development are among the most dramatic of the birth defects and have been quick to catch pubic attention ever since the thalidomide debacle of 1958 [1]. Indeed the tragedy of thalidomide provided the original impetus for the development of the modern system of pharmaceutical drug regulation and approval which has since been instituted worldwide.
Phocomelia (“flipper limb”) involves the absence of the upper two segments of the limb and the joining of the hand or foot to the proximal joint. Amelia (“no limb”) demonstrates complete limb absence. Two definitive studies have been published in recent years by the International Clearing House of Birth Defects and Surveillance Research [2, 3], which surveyed birth defects registries globally on phocomelia and amelia and found marked variation in the frequency of these anomalies across different nations globally. Curiously the only place that scored maximally in both reports was Melbourne in Australia. That in itself is a fascinating finding as Melbourne is a city with high numbers of immigrants from all over the world. The fact that their countries of origin, which are largely in Europe, have much lower rates of these anomalies points strongly to a local factor as an environmental teratogen. For several decades the city has been surrounded by cannabis farms [4–6].
Phocomelia was first noted to occur in experimental studies of white rabbits and hamsters following prenatal cannabis exposure [7–9]. Phocomelia [odds ratio (OR) = 21.90, 95% confidence interval (CI) 4.45–65.63], syndactyly (OR 24.33, 10.40–48.63) and polydactyly (OR 8.87, 4.82–14.87) were also identified to be elevated in an Hawaiian series looking at substance associations of prenatal drug exposures [10] with cannabis use identified as the leading association. Limb anomalies and leg reductions were also recently reported in a large national series from the USA [11, 12]. More recently elevated rates of limb reduction anomalies were also noted to have occurred in Europe [13].
In recent years, great strides have been made in the mechanistic understanding of both cannabis genotoxicity and limb development. Limb development has been shown to be controlled by an ordered and highly orchestrated sequence of events whereby proximal–distal, anterior–posterior and left–right gradients of embryonic tissue morphogens control limb bud development and outgrowth [14]. Interference with these delicate and finely balanced systems whether directly or epigenetically can clearly perturb, disrupt or completely truncate limb development either in part or in its totality depending on the timing of exposure.
CAs form one major clinical expression of cannabis genotoxicity. Genotoxic actions may also be expressed as elevations in rates of cancer or of cellular or organismal aging. Increased rates of cancer [15–20] and of aging [21] have all been described following cannabis exposure. For this reason, the present discussion is only part of the whole corpus of information relating to cannabis genotoxicity and it is important to bear the totality of the information in the field in mind.
A variety of cellular mechanisms have been described which might account for cannabis-induced genotoxicity including abnormalities of mitotic and meiotic cell divisions [22–26], sperm morphology and development [25], oocyte cell division [27], many abnormalities in the male and female reproductive tracts [28–30], single- and double-stranded DNA breaks [22, 24, 25, 31], chromosomal end-to-end fusions including ring and chain formation and the formation of double minute circles and micronuclei [25, 26], which are themselves engines for chromothriptic events [32] and oxidation of DNA bases [31]. The epigenome is also altered both in terms of DNA methylation [33–38] and histone synthesis and post-translational modification [39–44], and both DNA methylation changes [45, 46] and histone modifications [47] have been shown to be heritable via sperm. Cannabis inhibits telomerase [45, 48].
In terms of environmental contamination with known teratogens, it is a matter of record that both Europe and the USA have experienced a triple convergence of rising cannabis use, rising intensity of daily use and rising potency for Δ9-tetrahydrocannabinol (THC) concentration in both cannabis herb and resin in recent years in many of their regions [49–56]. This triple convergence driving cannabinoid exposure upward is exacerbated by the presumptive entry of cannabinoid and hemp products into the food chain of parts of Europe and the USA [49, 50] and by the well-documented exponential genotoxic effect of cannabinoids at higher doses which has been amply demonstrated in both laboratory [42, 57–62] and epidemiological studies [11–13, 15, 18–20, 63].
There is therefore a real concern that some communities are propelled into higher cannabinoid-use zones; this high-dose genotoxic exponentiation could be experienced as a “switch-like” mechanism for CAs as was recently reported from north eastern France where limb reduction defects at 60 times the historical controls were reported [64–66]. It is noted that the 60-fold increased incidence lies within the confidence interval (CI) reported from Hawaii [10]. Large cannabis crops are grown in the area [67]. It has been observed at the same time that French cows are also being born without limbs [64–66], thereby pointing directly to the food chain as the source of the environmental teratogen. Similar reports have also come from Germany [68] as cannabis use there rises [50]; however, they are not rising in nearby Switzerland where cannabis products are not permitted in the food chain.
It is important to note that close correlations between state and national levels of cannabis use and parental and maternal cannabis use have been demonstrated by multiple studies [69–76].
The present study builds on results recently reported in a major survey of European CAs [13] and explores these findings in depth. The hypotheses were generated prior to the commencement of the analysis. For all of these reasons, the present report investigated the substance and cannabis relationships of limb CAs in Europe in a formal time-lagged multivariable regression framework, in a causal inferential analytical paradigm and with regard to their native space–time context.
Methods
Data
Data on all available congenital anomaly rates (CARs) were downloaded by each individual year for each of 14 nations from the European Network of Population-Based Registries for the Epidemiological Surveillance of Congenital Anomalies (EUROCAT) website [77] and analyzed. EUROCAT total CAR includes anomaly rates among live births, stillbirths and cases where early termination for anomaly was practiced all combined together so that it represents a total overall picture across all classes of births. The nations selected were chosen on the basis of the availability of their CA data for most of the years 2010–2019. National tobacco (percentage daily tobacco use prevalence) and alcohol (liters of pure alcohol consumed per capita annually) use data were downloaded from the World Health Organization [78]. Drug use data for cannabis, amphetamines and cocaine was taken from the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) [79]. Last month cannabis use data were also supplemented by data on the THC content of cannabis herb and resin provided in recent published reports [50]. Data on daily cannabis use were also available from EMCDDA and were collated in recent reports [50]. Median household income data (in USD) were taken from the World Bank [80].
Anomaly definitions
The anomalies considered were: talipes (club foot), hip dysplasia, limb reduction anomalies (including amelia, hemimelia, phocomelia and their variants), polydactyly (extra fingers or toes) and syndactyly (fingers or toes fused together) and limb anomalies overall as a total group.
National assignment
Nations were categorized as being either high and rising daily cannabis use or low and/or falling daily cannabis use based on a recent European epidemiological study (see Supplementary Fig. 4) [50]. Thus, Belgium, Croatia, France, Germany, Italy, Netherlands, Norway, Portugal and Spain were categorized as nations experiencing increasing daily use, while Bulgaria, Finland, Hungary, Poland and Sweden were nations which were experiencing low or falling levels of daily cannabis use.
Derived data
The availability of several metrics of cannabis use, exposure and consumption made it possible to calculate various derived metrics. Hence, last month cannabis use prevalence data was multiplied by the THC content of cannabis herb and resin to derive compound metrics. These metrics were also multiplied by imputed daily cannabis use prevalence rates to derive further compound metrics for both cannabis herb and resin.
Data imputation
Missing data were completed by linear interpolation. This was particularly the case for daily cannabis use. Fifty-nine data points on daily cannabis use from EMCDDA were available for these 14 nations across this period. Linear interpolation expanded this dataset to 129 data points (further details provided in the Results section). Data on cannabis resin THC concentration were not available for Sweden. However, it was noted that the resin to herb THC concentration was almost constant in nearby Norway at 17.7 so this ratio was applied to the Swedish cannabis herb THC concentration data to derive estimates of Swedish cannabis resin THC concentration. Similarly, data for the cannabis resin THC concentration in Poland were not available. The resin to herb THC concentration ratio of nearby Germany was used to estimate the resin THC content in Poland from the known Polish herb THC concentrations. Since geospatial analytical techniques do not tolerate missing data the dataset was completed by the last observation carried forward or backward for Croatia in 2018 and 2019 and Netherlands in 2010. It was not appropriate to use multiple imputation methods for this dataset as multiple imputations cannot be applied in panel or spatial multivariable regression techniques.
Statistics
Data were processed in R Studio version 1.4.1717 based on R version 4.1.1 from the Comprehensive R Archive Network and the R Foundation for Statistical Computing [81]. The analysis was conducted in December 2021. Data were manipulated using dplyr from the tidyverse [82]. Data were log transformed where appropriate to improve compliance with normality assumptions based on the results of the Shapiro–Wilks test. Graphs were drawn in ggplot2 from tidyverse. Maps were drawn using ggplot2, sf (simple features) [83] and both custom color palettes and palettes taken from the viridis and viridisLite packages [84].
Bivariate maps were drawn with package colorplaner [85]. All illustrations are original and have not been published previously. Linear regression was conducted in Base R. Mixed effects regression was performed with the nation as the random effect using package nlme [86]. In all multivariable models, model reduction was by the classical technique of serial deletion of the least significant term to yield a final reduced model which is the model presented. Multiple linear models were processed in a single pass using combined techniques from R packages purrr and broom [82, 87, 88]. Different forms of regression were used to show that the observed effects were independent of the analytical method employed. The overall effect of covariates in multivariable models may be quantified as the marginal effect. In this case the overall marginal effect was calculated using the R package margins [89].
Categorical analyses
Cannabis exposure metrics can be categorized into exposure quintiles. This has been done from the lowest to highest quintiles numbering 1 to 5. These data can be paired with the relevant nations so that numbers exposure in each quintile can be computed. From the exposed numbers and their rates prevalence ratios (which are like relative rates in cohort studies), attributable fractions in the exposed and population attributable risks (also known as attributable fraction in the population) can be calculated. Preliminary analyses showed that, as might be expected, comparison of the highest and lowest quintiles gave the best separation between groups.
Covariate selection
The presence of multiple different metrics for cannabis consumption and exposure created a problem for analysis as it was not clear which was the most appropriate metric to employ for any particular model. Indiscriminate use of excessive covariates in a multivariable model would unnecessarily consume degrees of freedom and thereby restrict ability to assess interactions. This issue was formally addressed by the use of random forest regression using the R package ranger [90] with variable importance being formally assessed via the R package vip (variable importance plot) [91]. The most predictive covariates from this process were entered into the regression modeling equations. The tables from this analysis are presented in the Results section.
Panel and geospatial analysis
Panel analysis was conducted using R package plm [92] across both space and time simultaneously using the “twoways” effect. Panel analysis was conducted because it lends itself to lagged temporal analyses and can also be inverse probability weighted. Inverse probability weighting is not possible with extant geospatial modeling techniques at this time. The spatial weights matrix was calculated using the edge and corner “queen” relationships using R package spdep (spatial dependency) [93]. Geospatial modeling was conducted using the spatial panel random effects maximum likelihood (spreml) function from the package spml which allows detailed modeling and correction of model error structures [94, 95]. Such models may produce four model coefficients of interest which are useful in determining the most appropriate error structure for the model. These coefficients are phi the random error effect, psi the serial correlation effect, rho the spatial coefficient and theta the spatial autocorrelation coefficient. In each case the most appropriate error structure was chosen for each spatial model, generally taking care to preserve the model error specification across related models. The appropriate error structure was determined by the backward methods from the full general model to the most specific model as has been described [96]. Both panel and geospatial models were temporally lagged as indicated by 1–2 years.
Causal inference
The formal tools of causal inference were used in this analysis. Inverse probability weighting (ipw) is the technique of choice to convert a purely observational study into a pseudo-randomized study from which it is appropriate to make causal inferences [97]. All multivariable panel models presented herein were inverse probability weighted. Inverse probability weighting was performed using the R package ipw. Similarly E-values (expected values) quantify the correlation required of some hypothetical unmeasured confounder covariate with both the exposure of concern and the outcome of interest in order to explain away some apparently causal relationship [98–100]. It therefore provides a quantitative measure of the robustness of the model to extraneous covariates which have not been accounted for within the measured parameters. E-values represent a transformation of the relative risk and can be calculated in several ways such as from the raw numbers in the exposed and control groups or from the output of regression models where the model standard deviation is known. E-values have a CI associated with them and the 95% lower bound of this CI is reported herein. E-value estimates greater than 1.25 are said to indicate causality [101] with E-values greater than nine being described as high [102]. E-values were calculated from the R package EValue [103]. Both inverse probability weighting and E-values are foundational and pivotal techniques used in formal causal inferential methods in order to allow causal relationships to be assessed from real-world observational studies.
Data availability
Raw datasets including 3800 lines of computation code in R has been made freely available through the Mendeley data repository at the following URLs: 10.17632/vsmmmkncsd.1 and 10.17632/vd6mt5r5jm.1.
Ethics
Ethical approval for this study was provided from the Human Research Ethics Committee of the University of Western Australia number RA/4/20/4724 on 24 September 2021.
Results
The overall plan of the presentation of the results of this analysis is to first present univariate data, then bivariate continuous and categorical data, then move progressively to multivariable regression by mixed effects, panel and geospatial techniques. Finally, a detailed consideration of E-values from both specific models and overall analyses and their implications for causal inference is presented.
Seven-hundred and two CARs for limb-related birth defects are obtained for the 14 nations listed in Supplementary Table ST1. As shown in this table, 107–122 CARs were obtained in each of the six groups of limb defects in general, limb reduction defects, talipes (club foot), hip dislocation, polydactyly and syndactyly. The terminology here is somewhat confusing in that both whole classes of anomalies are referred to as limb anomalies, but also one specific anomaly in particular is also referred to as a limb anomaly. To deconfuse this issue in the present report when reference is made to the class of limb anomalies the group of six anomalies will be referred to as the class limb anomalies. The term limb anomaly will be reserved for the limb anomalies as a specific sub-group of the class of anomalies.
Drug exposure rates including various compound metrics of cannabis exposure are as shown in ST1 along with national median household income.
As discussed in the Methods section the rates of daily cannabis use were missing in many years by jurisdiction. The 59 raw data points are shown in Supplementary Table ST2 and further 70 data points were completed by linear interpolation (a total of 129 data points altogether) are shown in Supplementary Table ST3.
As noted above Belgium, Croatia, France, Germany, Italy, Netherlands, Norway, Portugal and Spain were categorized as nations experiencing increasing daily cannabis use, while Bulgaria, Finland, Hungary, Poland and Sweden were nations which were experiencing low or falling levels of daily cannabis use. Figure 1 shows the trend lines for the rates of these six anomalies according to the five substances tobacco, alcohol, cannabis resin THC content, amphetamine and cocaine exposure. While most of these substances demonstrate a weak or overtly negative relationship with CARs, it is apparent that the relationship between cannabis resin and CAR is relatively strong and positive (noted by the rising slope of the regression line in each case).
Figure 2 plays a similar role showing the CARs against the various cannabis exposure metrics. Overall, the relationship for the resin THC content and the daily cannabis use interpolated seem to be most strongly related metrics.
As noted in the Methods section the 14 nations were divided into countries where daily cannabis use is high and or rising and those where it is low and or decreasing [50]. In Bulgaria, Finland, Hungary, Poland and Sweden the rate of daily cannabis use was declining while it was high or rising in the other countries. Comparing the CARs between the two sets of nations categorized in this way produced the appearances shown in Fig. 3 where countries with rising rates of cannabis use appear to have higher rate of hip dysplasia and limb anomalies than the other group of countries. At linear regression the rate of limb anomalies in the nations with increasing cannabis use was significantly higher than those without (β-est. = 0.3823, t = 3.58, P = 0.0005; model F = 12.8, df = 1, 105, P = 0.0005) and similarly for hip dysplasia (β-est. = 2.054, t = 9.7911, P = 1.81 × 10−16; model F = 95.86, df = 1, 105, P = 1.81 × 10−16). At mixed effects regression across all six anomalies using the anomaly as the random effect nations with increasing daily use had higher overall limb CARs (β-est. = 1.502 × 10−4, t = 4.745, P = 2.53x10−6; model AIC = 1.61 × 103, LogLik. = −801.40).
Figure 4 shows the progression of limb anomaly CARs across Europe over time. High rates in France and low rates in some countries are immediately apparent. The pattern of limb reduction anomalies is shown in Fig. 5. Supplementary Figures SF1–SF3 show a similar role for polydactyly, syndactyly and hip dysplasia. Polydactyly rates are relatively high in France. Syndactyly rates are high in Finland and Hungary but then their data are unavailable. Earlier in the time series, hip dysplasia rates are high in Finland but then France and Netherlands supersede the Finns.
Supplementary Figure SF4 shows the rate of the compound index which is the product of last month cannabis use, cannabis resin THC concentration and interpolated daily use. This index seems to be highest in Spain and France where the index has risen across the decade. Figure 6 shows the rate of daily interpolated use across Europe. Again rises in Netherlands, France and Spain are noteworthy.
Supplementary Figure SF5 is a bivariate map graph of the bivariate relationship between limb anomalies and cannabis herb THC concentration. Countries which are shaded green have both low cannabis herb THC concentration and low limb anomaly rates. Countries which are shaded bright pink or purple have high rates of both covariates. Hence, across time France changes from red to bright pink and Germany stays purple indicating convergence of high rates of cannabis herb THC concentration and limb anomaly rates. The meaning of other colors is as shown in the colorplane key.
Figure 7 performs a similar role for the relationship between cannabis resin THC concentration and limb anomaly rates. In this study, France, Germany, Netherlands and Belgium are noted to turn from reddish to pink and purple across the time series.
Figure 8 performs a similar role for limb reduction anomalies as a function of cannabis resin THC concentration. Here France, Germany, Belgium, Spain and Bulgaria turn from reddish or green to pink or purple, showing a temporal coincidence between rising rates of both variables. Polydactyly (Fig. 9) and syndactyly (Fig. 10) show a broadly similar pattern.
In Fig. 11, France (including Corsica) is noted to turn pink as the incidence of hip dysplasia and the compound index last month cannabis use × resin THC concentration × daily use interpolated increase together. When hip dysplasia is considered against the compound index of last month cannabis use × resin THC concentration, France, Portugal and Netherlands are noted to turn pink and purple (Supplementary Fig. SF6). When hip dysplasia is considered against daily cannabis use interpolated, France turns somewhat pink and Netherlands turns purple (Supplementary Fig. SF7).
Supplementary Table ST4 sets out the formal assessment of these bivariate relationships by linear regression where all models and all substances are considered simultaneously through a purrr–broom analytical pathway. As well as the usual regression parameters, the E-value estimate and the lower 95% bound of its CI are also listed in the columns at the right side of the table. Models are listed in descending order of minimum E-value (mEV). It is noted that the most significant models which appear at the top of the table are dominated by various metrics of cannabis use and exposure.
Table 1 extracts the most significant models from the preceding table in descending order of mEV which are significantly elevated. Here 19 of the 22 terms listed feature various metrics of cannabis exposure. Eight of the 22 are above nine and so would be considered high [102], and 21/22 exceed 1.25 which is the generally quoted threshold for causality [101].
| Anomaly | Substance | Mean anomaly rate | Estimate | Std. error | Sigma | t_statistic | P_Value | E-value Estimate | E-value lower bound |
|---|---|---|---|---|---|---|---|---|---|
| Hip dysplasia | Daily.Interpol. | 5.7972 | 39.9028 | 12.1854 | 1.1146 | 3.2746 | 0.0014 | 2.82E + 14 | 9.96E + 05 |
| Limb | Daily.Interpol. | 36.3134 | 13.9285 | 4.8044 | 0.4394 | 2.8991 | 0.0046 | 6.72E + 12 | 2.38E + 04 |
| Limb | Herb | 36.3134 | 6.3467 | 1.2321 | 0.4080 | 5.1512 | 1.22E − 06 | 2.81E + 06 | 1.30E + 04 |
| Limb | LMCannabis_Herb | 36.3134 | 5.9187 | 1.6667 | 0.4315 | 3.5512 | 5.76E − 04 | 5.27E + 05 | 543.88 |
| Club foot | Herb | 10.5762 | 6.5507 | 2.1341 | 0.8298 | 3.0695 | 0.0027 | 2.64E + 03 | 26.58 |
| Club foot | Daily.Interpol. | 10.5762 | 20.4087 | 9.1947 | 0.8630 | 2.2196 | 0.0284 | 4.43E + 09 | 25.24 |
| Limb | Resin | 36.3134 | 1.8357 | 0.5015 | 0.4318 | 3.6603 | 4.14E − 04 | 95.28 | 11.59 |
| Polydactyly | Herb | 10.5235 | 3.6119 | 1.3516 | 0.5256 | 2.6723 | 0.0086 | 1.04E + 03 | 10.17 |
| Syndactyly | Resin | 4.7039 | 2.1677 | 0.7134 | 0.6413 | 3.0383 | 0.0030 | 42.82 | 5.43 |
| Limb | LMCannabis_Resin | 36.3134 | 1.2041 | 0.3842 | 0.4391 | 3.1344 | 0.0023 | 23.75 | 4.55 |
| Polydactyly | LMCannabis_Herb | 10.5235 | 3.9646 | 1.8495 | 0.5309 | 2.1436 | 0.0341 | 1.79E + 03 | 3.03 |
| Hip dysplasia | LM_Cannabis | 5.7972 | 11.1612 | 5.3423 | 1.1465 | 2.0892 | 0.0391 | 1.41E + 04 | 2.91 |
| Hip dysplasia | Cocaine | 5.7972 | 0.7591 | 0.1383 | 1.0314 | 5.4887 | 2.83E − 07 | 3.32 | 2.45 |
| Limb | Cocaine | 36.3134 | 0.2917 | 0.0542 | 0.4044 | 5.3797 | 4.56E − 07 | 3.27 | 2.41 |
| Polydactyly | LMCannabis_Resin | 10.5235 | 1.0689 | 0.4423 | 0.5296 | 2.4167 | 0.0174 | 12.03 | 2.19 |
| Club foot | LMCannabis_Herb | 10.5762 | 6.0400 | 2.9511 | 0.8471 | 2.0467 | 0.0429 | 1.31E + 03 | 2.00 |
| Limb | LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | 36.3134 | 1.1656 | 0.5302 | 0.4465 | 2.1982 | 0.0301 | 21.00 | 1.92 |
| Limb | LM.Cannabis_x_Resin.THC_x_Daily.Interpol. | 36.3134 | 0.5944 | 0.2401 | 0.4470 | 2.4752 | 0.0151 | 6.16 | 1.90 |
| Polydactyly | Resin | 10.5235 | 1.3047 | 0.5918 | 0.5320 | 2.2047 | 0.0297 | 18.12 | 1.89 |
| Club foot | Cocaine | 10.5762 | 0.3514 | 0.1024 | 0.8224 | 3.4315 | 8.24E − 04 | 2.31 | 1.65 |
| Club foot | LM.Cannabis_x_Resin.THC_x_Daily.Interpol. | 10.5762 | 1.0388 | 0.4812 | 0.9004 | 2.1589 | 0.0333 | 5.16 | 1.44 |
| Polydactyly | Daily.Interpol. | 10.5235 | 0.1369 | 0.0662 | 0.5316 | 2.0677 | 0.0408 | 1.84 | 1.13 |
Figure 12 presents a boxplot comparing the logarithm of the anomaly rate by the extreme highest and lowest quintiles of daily cannabis exposure interpolated. Non-overlapping notches on the boxes signify statistically significant differences between the two groups. For club foot, hip dysplasias, limb reductions and polydactyly, the highest quintile clearly experiences rates significantly above those of the lowest quintile. Figure 13 presents similar boxplots for the compound indices of daily cannabis use interpolated:last month cannabis use:resin THC concentrations and shows that the rates for club foot, hip dysplasia, limb anomalies and polydactyly are significant greater in the higher quintiles. Figure 14 presents boxplots for daily cannabis use interpolated:last month cannabis use:herb THC concentrations and finds that rates of club foot, hip dysplasias and limb anomalies are higher in the highest quintiles.
The numbers in the highest and lowest quintile groups together with their applicable prevalence ratios are indicated in Supplementary Table S5. These categorical comparisons are shown in tabular format in Table 2. This table also shows the attributable fraction in the exposed and the population attributable risk (attributable fraction in the population). Impressively low P-values for many anomaly groups especially for limb anomalies are observed.
| Calculated rates | Significance | E-values | ||||||
|---|---|---|---|---|---|---|---|---|
| Congenital anomaly | Covariate | Prevalence ratio (CI) | Attributable fraction in the exposed (CI) | Population attributable risk (CI) | Chi-squared | P-value | E-value—Point | E-value—Lower |
| Club foot—talipes equinovarus | Daily_Use | 2.0564 (1.9767, 2.1394) | 0.5137 (0.4941, 0.5326) | 0.4394 (0.4201, 0.458) | 1332.79 | 4.23E − 292 | 3.53 | 3.37 |
| Club foot—talipes equinovarus | pm.Resin.Daily | 1.8754 (1.8032, 1.9505) | 0.4668 (0.4454, 0.4873) | 0.3875 (0.3672, 0.4071) | 1018.44 | 8.81E − 224 | 3.16 | 3.01 |
| Club foot—talipes equinovarus | Herb_THC | 1.7624 (1.694, 1.8336) | 0.4326 (0.4097, 0.4546) | 0.3619 (0.3404, 0.3827) | 807.62 | 5.95E − 178 | 2.92 | 2.78 |
| Club foot—talipes equinovarus | pm.Herb.Daily | 1.8163 (1.7403, 1.8957) | 0.4494 (0.4254, 0.4725) | 0.3911 (0.368, 0.4133) | 770.93 | 5.63E − 170 | 3.03 | 2.88 |
| Club foot—talipes equinovarus | LM.Cann_Resin_THC | 1.5491 (1.5007, 1.5991) | 0.3545 (0.3336, 0.3747) | 0.2523 (0.2352, 0.269) | 741.04 | 1.78E − 163 | 2.47 | 2.37 |
| Club foot—talipes equinovarus | LM.Cann_Herb_THC | 1.482 (1.4226, 1.5439) | 0.3252 (0.2971, 0.3523) | 0.2785 (0.2528, 0.3033) | 360.02 | 1.39E − 80 | 2.33 | 2.20 |
| Club foot—talipes equinovarus | Resin_THC | 1.4277 (1.3726, 1.485) | 0.2996 (0.2715, 0.3266) | 0.1147 (0.1013, 0.128) | 317.85 | 2.13E − 71 | 2.21 | 2.09 |
| Club foot—talipes equinovarus | LM_Cannabis | 1.3922 (1.341, 1.4454) | 0.2817 (0.2543, 0.3082) | 0.235 (0.2108, 0.2586) | 302.29 | 5.22E − 68 | 2.13 | 2.02 |
| Hip dysplasia | Daily_Use | 12.0452 (10.8119, 13.4191) | 0.917 (0.9075, 0.9255) | 0.8931 (0.8813, 0.9038) | 3336.94 | 0.00E + 00 | 23.58 | 21.11 |
| Hip dysplasia | pm.Herb.Daily | 7.2802 (6.5544, 8.0864) | 0.8626 (0.8474, 0.8763) | 0.8367 (0.8192, 0.8525) | 1887.95 | 0.00E + 00 | 14.04 | 12.59 |
| Hip dysplasia | pm.Resin.Daily | 3.0124 (2.8457, 3.1889) | 0.668 (0.6486, 0.6864) | 0.5972 (0.5762, 0.6172) | 1593.92 | 0.00E + 00 | 5.47 | 5.14 |
| Hip dysplasia | Resin_THC | 2.4282 (2.3027, 2.5604) | 0.5882 (0.5657, 0.6094) | 0.3029 (0.2836, 0.3216) | 1146.86 | 1.08E − 251 | 4.29 | 4.03 |
| Hip dysplasia | LM_Cannabis | 2.4663 (2.3009, 2.6435) | 0.5945 (0.5654, 0.6217) | 0.5537 (0.5239, 0.5816) | 695.46 | 1.45E − 153 | 4.37 | 4.03 |
| Hip dysplasia | LM.Cann_Herb_THC | 2.4583 (2.2555, 2.6793) | 0.5932 (0.5566, 0.6268) | 0.5628 (0.5256, 0.5971) | 448.49 | 7.70E − 100 | 4.35 | 3.94 |
| Hip dysplasia | Herb_THC | 1.4499 (1.3648, 1.5402) | 0.3103 (0.2673, 0.3507) | 0.2625 (0.2238, 0.2993) | 146.78 | 4.38E − 34 | 2.26 | 2.07 |
| Hip dysplasia | LM.Cann_Resin_THC | 2.4249 (2.3135, 2.5416) | 0.5876 (0.5678, 0.6065) | 0.4832 (0.4628, 0.5027) | 1455.89 | 7.51E − 319 | 4.28 | 4.06 |
| Limb | LM_Cannabis | 1.7369 (1.6932, 1.7816) | 0.4242 (0.4094, 0.4387) | 0.3840 (0.3696, 0.3980) | 1856.66 | 0.00E + 00 | 2.87 | 2.78 |
| Limb | Herb_THC | 1.7137 (1.6716, 1.7568) | 0.4165 (0.4018, 0.4308) | 0.3609 (0.3470, 0.3745) | 1848.39 | 0.00E + 00 | 2.82 | 2.73 |
| Limb | Resin_THC | 1.5997 (1.5669, 1.6332) | 0.3749 (0.3618, 0.3877) | 0.1543 (0.1470, 0.1615) | 2014.71 | 0.00E + 00 | 2.58 | 2.51 |
| Limb | Daily_Use | 2.1803 (2.1337, 2.2279) | 0.5413 (0.5313, 0.5512) | 0.4718 (0.4618, 0.4817) | 5265.30 | 0.00E + 00 | 3.78 | 3.69 |
| Limb | LM.Cann_Herb_THC | 2.0594 (1.9927, 2.1283) | 0.5144 (0.4982, 0.5302) | 0.4833 (0.4670, 0.4990) | 1935.56 | 0.00E + 00 | 3.54 | 3.40 |
| Limb | LM.Cann_Resin_THC | 1.6726 (1.6428, 1.7028) | 0.4021 (0.3913, 0.4127) | 0.3061 (0.2966, 0.3156) | 3233.77 | 0.00E + 00 | 2.73 | 2.67 |
| Limb | pm.Herb.Daily | 2.2839 (2.2246, 2.3449) | 0.5622 (0.5505, 0.5735) | 0.5117 (0.4998, 0.5232) | 4007.14 | 0.00E + 00 | 4.00 | 3.88 |
| Limb | pm.Resin.Daily | 2.0808 (2.0362, 2.1265) | 0.5194 (0.5089, 0.5297) | 0.4433 (0.4329, 0.4536) | 4590.80 | 0.00E + 00 | 3.58 | 3.49 |
| Limb reductions | Daily_Use | 1.5446 (1.4684, 1.6249) | 0.3526 (0.3190, 0.3846) | 0.2878 (0.2577, 0.3166) | 287.62 | 8.20E − 65 | 2.46 | 2.30 |
| Limb reductions | pm.Resin.Daily | 1.5247 (1.4449, 1.6089) | 0.3441 (0.3079, 0.3785) | 0.2749 (0.2431, 0.3054) | 240.02 | 1.94E − 54 | 2.42 | 2.25 |
| Limb reductions | pm.Herb.Daily | 1.4774 (1.3965, 1.5631) | 0.3231 (0.2839, 0.3602) | 0.2731 (0.2376, 0.3069) | 186.66 | 8.51E − 43 | 2.32 | 2.14 |
| Limb reductions | LM.Cann_Resin_THC | 1.3277 (1.2702, 1.3878) | 0.2468 (0.2127, 0.2794) | 0.1676 (0.1422, 0.1922) | 158.70 | 1.09E − 36 | 1.99 | 1.86 |
| Limb reductions | LM_Cannabis | 1.3606 (1.2892, 1.4360) | 0.2650 (0.2243, 0.3036) | 0.2203 (0.1845, 0.2544) | 126.28 | 1.34E − 29 | 2.06 | 1.90 |
| Limb reductions | Herb_THC | 1.3302 (1.2608, 1.4034) | 0.2482 (0.2069, 0.2875) | 0.1972 (0.1623, 0.2307) | 109.66 | 5.80E − 26 | 1.99 | 1.83 |
| Limb reductions | LM.Cann_Herb_THC | 1.2620 (1.1921, 1.3360) | 0.2076 (0.1612, 0.2515) | 0.1734 (0.1331, 0.2119) | 64.39 | 5.10E − 16 | 1.84 | 1.67 |
| Limb reductions | Resin_THC | 1.1734 (1.1109, 1.2395) | 0.1478 (0.0998, 0.1932) | 0.0499 (0.0322, 0.0674) | 32.81 | 5.08E − 09 | 1.62 | 1.46 |
| Polydactyly | pm.Resin.Daily | 2.0453 (1.9666, 2.1271) | 0.5111 (0.4915, 0.5299) | 0.4303 (0.4112, 0.4488) | 1334.92 | 1.46E − 292 | 3.51 | 3.35 |
| Polydactyly | LM.Cann_Resin_THC | 1.6572 (1.6067, 1.7092) | 0.3966 (0.3776, 0.4149) | 0.2876 (0.2715, 0.3034) | 1047.40 | 4.46E − 230 | 2.70 | 2.59 |
| Polydactyly | Daily_Use | 1.7222 (1.6613, 1.7853) | 0.4193 (0.398, 0.4399) | 0.3489 (0.329, 0.3681) | 897.24 | 1.95E − 197 | 2.84 | 2.71 |
| Polydactyly | Resin_THC | 1.6629 (1.6021, 1.7259) | 0.3986 (0.3758, 0.4206) | 0.1673 (0.1542, 0.1802) | 732.61 | 1.21E − 161 | 2.71 | 2.58 |
| Polydactyly | pm.Herb.Daily | 1.6546 (1.5909, 1.7207) | 0.3956 (0.3714, 0.4189) | 0.3399 (0.3173, 0.3618) | 646.99 | 5.03E − 143 | 2.70 | 2.56 |
| Polydactyly | LM_Cannabis | 1.5813 (1.5227, 1.6422) | 0.3676 (0.3433, 0.3911) | 0.3129 (0.2904, 0.3346) | 575.34 | 1.94E − 127 | 2.54 | 2.41 |
| Polydactyly | Herb_THC | 1.611 (1.5471, 1.6775) | 0.3793 (0.3536, 0.4039) | 0.3125 (0.2892, 0.335) | 544.05 | 1.24E − 120 | 2.60 | 2.47 |
| Polydactyly | LM.Cann_Herb_THC | 1.5969 (1.5335, 1.6629) | 0.3738 (0.3479, 0.3987) | 0.3234 (0.2993, 0.3467) | 522.71 | 5.44E − 116 | 2.57 | 2.44 |
| Syndactyly | pm.Resin.Daily | 1.2841 (1.2082, 1.3647) | 0.2212 (0.1724, 0.2672) | 0.1703 (0.1305, 0.2083) | 65.15 | 3.47E − 16 | 1.89 | 1.71 |
| Syndactyly | Resin_THC | 1.1986 (1.1223, 1.2801) | 0.1657 (0.109, 0.2188) | 0.0568 (0.0353, 0.0778) | 29.20 | 3.27E − 08 | 1.69 | 1.49 |
| Syndactyly | Herb_THC | 1.0751 (1.0142, 1.1398) | 0.0699 (0.014, 0.1226) | 0.0529 (0.0101, 0.0939) | 5.92 | 0.0075 | 1.36 | 1.13 |
Supplementary Table ST6 presents a similar tabulation for an additive panel model which has been inverse probability weighted. Here 15 of the 48 terms listed demonstrate elevated mEVs. The first four terms in this table are all occupied by cannabis resin THC concentration, consistent with the data appearances presented in Fig. 2. This table also includes a calculation of the average marginal effect contrasting the lowest level exposure to the highest level exposure (minimum − maximum comparison) in both absolute terms on the scale of the response variable, CAR and as a relative percentage of the CAR.
In order to see if the results would be robust to the use of another modeling technique, we also employed mixed effects regression with inverse probability weighting. Table 3 presents the results of serialized mixed effects models in purrr and lists the models in order of the CA and descending mEV. Limb anomalies, limb reduction anomalies, polydactyly and syndactyly are all significantly related to cannabis resin THC concentration and talipes is related to cannabis herb THC concentration. Supplementary Table ST7 lists these terms in order of the covariate assessed. Annotations for the THC concentration of both cannabis resin and herb again feature.
| Anomaly | Term | Mean rate | Estimate | Std. error | Sigma | Adj. R squared | t_statistic | P_Value | AME—Min. to Max. | E-value estimate | E-value lower bound |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Club foot | Herb | 10.647 | 6.514 | 2.755 | 0.684 | 0.378 | 2.364 | 0.0202 | 0.931 | 11 588.01 | 8.39 |
| Club foot | Alcohol | 10.647 | 0.097 | 0.049 | 0.684 | 0.378 | 1.993 | 0.0492 | 0.606 | 1.53 | 1.05 |
| Club foot | Income | 10.647 | 1.63E − 05 | 7.97E − 06 | 0.684 | 0.378 | 2.048 | 0.0434 | 1.225 | 1.00 | 1.00 |
| Hip dysplasia | Cocaine | 6.349 | 0.824 | 0.331 | 1.167 | 0.268 | 2.489 | 0.0147 | 2.048 | 3.21 | 1.56 |
| Limb | Resin | 37.746 | 2.128 | 0.628 | 0.478 | 0.484 | 3.387 | 0.0011 | 1.097 | 114.38 | 10.55 |
| Limb | Cocaine | 37.746 | 0.541 | 0.136 | 0.478 | 0.484 | 3.990 | 1.36E − 04 | 1.345 | 5.05 | 2.77 |
| Limb | Tobacco | 37.746 | 0.058 | 0.018 | 0.478 | 0.484 | 3.144 | 0.0023 | 1.438 | 1.48 | 1.25 |
| Limb | Alcohol | 37.746 | 0.071 | 0.034 | 0.478 | 0.484 | 2.076 | 0.0408 | 0.441 | 1.55 | 1.10 |
| Limb | Income | 37.746 | 2.12E − 05 | 5.57E − 06 | 0.478 | 0.484 | 3.799 | 2.67E − 04 | 1.588 | 1.01 | 1.00 |
| Limb reductions | Resin | 5.389 | 3.822 | 0.822 | 0.626 | 0.080 | 4.651 | 1.12E − 05 | 1.970 | 517.93 | 49.56 |
| Limb reductions | Cocaine | 5.389 | 0.595 | 0.174 | 0.626 | 0.080 | 3.429 | 9.12E − 04 | 1.479 | 4.18 | 2.26 |
| Limb reductions | Alcohol | 5.389 | 0.153 | 0.045 | 0.626 | 0.080 | 3.424 | 9.26E − 04 | 0.952 | 1.81 | 1.43 |
| Limb reductions | Tobacco | 5.389 | 0.098 | 0.024 | 0.626 | 0.080 | 4.039 | 1.12E − 04 | 2.418 | 1.57 | 1.36 |
| Limb reductions | Income | 5.389 | 2.74E − 05 | 7.29E − 06 | 0.626 | 0.080 | 3.764 | 2.96E − 04 | 2.059 | 1.01 | 1.00 |
| Polydactyly | Cocaine | 10.350 | 1.043 | 0.163 | 0.588 | 0.240 | 6.401 | 6.59E − 09 | 2.593 | 9.53 | 5.59 |
| Polydactyly | Resin | 10.350 | 2.197 | 0.772 | 0.588 | 0.240 | 2.847 | 0.0055 | 1.133 | 59.55 | 5.25 |
| Polydactyly | Tobacco | 10.350 | 0.061 | 0.023 | 0.588 | 0.240 | 2.685 | 0.0086 | 1.509 | 1.43 | 1.19 |
| Polydactyly | Alcohol | 10.350 | 0.092 | 0.042 | 0.588 | 0.240 | 2.206 | 0.0299 | 0.576 | 1.58 | 1.15 |
| Syndactyly | Resin | 4.230 | 2.160 | 0.862 | 0.656 | 0.033 | 2.507 | 0.0140 | 1.114 | 39.48 | 3.27 |
| Syndactyly | Cocaine | 4.230 | 0.537 | 0.182 | 0.656 | 0.033 | 2.952 | 0.0040 | 1.335 | 3.63 | 1.89 |
Some of these metrics from mixed effects models are then summarized in Table 4 again in descending order of mEV. To make it easier to understand this table the data are also presented graphically in Figs 15 and 16. Figure 15 lists the number of anomalies associated with each exposure variable, the sum and mean of the negative P-value exponents in the top row, and the total, mean and median of the percentage marginal effect from the minimum to the maximum range of the covariate along the bottom row. In the top row of graphs resin THC concentration occupies the fifth or third position along the x-axis. It also occupies a largely mid-position rank in the lower set of graphs. Figure 16 presents an analysis of the sum, mean and median of the E-value estimates (top row) and the minimum E-values (bottom row). In five of the six cases cannabis herb and resin THC concentration occupy the top two positions along the x-axis.
| Term | Number anomalies | Sum P exponent | Mean P exponent | Median P exponent | Sum E-value Est. | Mean E-value Est. | Median E-value Est. | Sum E-value lower | Mean E-value lower | Median E-value lower | Total % increase AME Min − Max | Mean % increase AME Min − Max | Median % increase AME Min − Max |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Herb_THC | 1 | 1 | 1 | 1 | 11 588 | 11 588.01 | 11 588 | 8.39 | 8.39 | 8.39 | 10.07 | 10.07 | 10.07 |
| Resin_THC | 4 | 9 | 2.25 | 2 | 731.35 | 182.84 | 86.96 | 68.63 | 17.16 | 7.90 | 114.29 | 28.57 | 22.77 |
| Cocaine | 5 | 18 | 3.6 | 3 | 25.61 | 5.12 | 4.18 | 14.07 | 2.81 | 2.26 | 209.43 | 41.89 | 41.82 |
| Tobacco | 3 | 7 | 2.33 | 2 | 4.48 | 1.49 | 1.48 | 3.8 | 1.27 | 1.25 | 129.31 | 43.1 | 20.78 |
| Alcohol | 4 | 6 | 1.5 | 1 | 6.47 | 1.62 | 1.56 | 4.73 | 1.18 | 1.12 | 33.58 | 8.4 | 5.96 |
| Income | 3 | 7 | 2.33 | 3 | 3.02 | 1.01 | 1.01 | 3.01 | 1 | 1 | 92.37 | 30.79 | 14.61 |
The question then arises as to the most appropriate set of covariates to use in multiple regression models for each anomaly. This was determined by the use of random forest regression in conjunction with formal assessment of variable importance using the R packages ranger and vip together. Tables which rank the 13 covariates for each anomaly of interest are shown in Supplementary Tables ST8—ST12.
For the reasons explained in the Discussion section, we decided to focus our more detailed multivariable study on the five CAs of limb anomalies, limb reduction anomalies, hip dysplasia, polydactyly and syndactyly.
Supplementary Table ST13 sets out three inverse probability weight panel regression models for the congenital limb anomalies for additive, interactive and interactive models lagged by 2 years respectively. Details of model specifications are provided in each model in the table. It is noted that in each model terms including cannabis metrics are significant and positive.
Supplementary Table ST14 sets out final inverse probability weighted panel models for limb reduction anomalies. In this series of multivariable models, terms including cannabis are noted to be significant and positive in additive, interactive and interactive models lagged by 1 year, but this signal is lost after 2 years of temporal lag.
It is possible that, like the relationship between tobacco and lung cancer, the relationship between various cannabis metrics and limb reduction may be a polynomial function of exposure to environmental intoxicants. Indeed the tobacco–lung cancer relationship is known to be a fourth-order relationship [104]. For this reason various polynomial relationships were tested with the results of final models shown in Supplementary Table ST15. Quadratic, cubic and quartic models in exposure to various cannabis metrics were studied. In the second third and fourth models the significance levels were increased over the interactive model presented in ST14.
Additive, interactive and lagged inverse probability weighted panel models for hip dysplasia are presented in Supplementary Table ST16. In both the second and third models terms including cannabis metrics are positive and significant.
A similar pattern is continued in Supplementary Tables ST17 and 18 for polydactyly and syndactyly where in all models tested cannabis metrics persist in final models, have positive coefficients and are highly statistically significant.
A spatial weights matrix was constructed for the nations having complete or almost complete CA datasets. The last observation carried forward technique was used to complete missing datasets for Netherlands in 2010 and Croatia in 2018 and 2019 as described in the Methods section. The spatial weights matrix was constructed using the dedicated R package spdep. The geospatial links which form the basis of the sparse spatial weights matrix in their native, edited and final forms are illustrated by the maps shown in Supplementary Fig. SF8.
Table 5 shows the final reduced geospatial models for limb anomalies (the CA) for additive, interactive and lagged models. Once again in each case terms for various cannabis metrics are positive and (usually) highly statistically significant.
| Parameter values | Model parameters | ||||
|---|---|---|---|---|---|
| Parameter | Estimate (CI) | P-value | Parameter | Value | Significance |
| Additive | |||||
| Rate ∼ Tobacco + Alcohol + Herb + Resin + Daily.Interpol. + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.04 (0.01, 0.06) | 0.0105 | phi | 2.6051 | 0.0761 |
| Herb | 1.83 (0.12, 3.54) | 0.0364 | rho | 0.6275 | 1.04E − 09 |
| Resin | 1.06 (0.46, 1.66) | 0.0006 | lambda | −0.5205 | 1.70E − 05 |
| Income | 0 (0, 0) | 1.51E − 06 | |||
| Interactive | |||||
| Rate ∼ Tobacco + Alcohol + Herb * LM.Cannabis_x_Resin.THC_x_Daily.Interpol. + Daily.Interpol. + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.05 (0.04, 0.07) | 1.38E − 11 | rho | −0.6317 | 1.66E − 09 |
| Herb | 1.82 (0.56, 3.08) | 0.0047 | lambda | 0.5449 | 5.85E − 09 |
| Cocaine | 0.15 (0.07, 0.23) | 0.0001 | |||
| Income | 0 (0, 0) | <2.2E − 16 | |||
| 2 Lags | |||||
| Rate ∼ Tobacco + Alcohol + Herb * LM.Cannabis_x_Resin.THC_x_Daily.Interpol. + Daily.Interpol. + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.05 (0.04, 0.07) | 5.03E − 10 | rho | −0.6228 | 2.77E − 08 |
| Herb | 2.09 (0.58, 3.6) | 0.0066 | lambda | 0.3986 | 0.000516 |
| Cocaine | 0.16 (0.07, 0.25) | 0.0007 | |||
| Income | 0 (0, 0) | 2.68E − 14 | |||
Geospatial models for limb reduction are presented in Table 6. In this series of models there is a positive signal at 1 year of temporal lag for the compound metric last month cannabis use × resin THC concentration but not at other time points. This may imply that the signal is positive at this time point but disappears thereafter.
| Parameter values | Model parameters | ||||
|---|---|---|---|---|---|
| Parameter | Estimate (CI) | P-value | Parameter | Value | Significance |
| Additive | |||||
| Rate ∼ Tobacco + Alcohol + Resin + Herb + LM.Cannabis_x_Resin.THC + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.07 (0.04, 0.09) | 1.11E − 07 | rho | −0.6489 | 1.96E − 09 |
| Income | 0 (0, 0) | 8.56E − 08 | lambda | 0.6118 | 2.38E − 09 |
| Interactive | |||||
| Rate ∼ Tobacco + Alcohol + Resin * Herb + LM.Cannabis_x_Resin.THC + + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.07 (0.04, 0.09) | 1.11E − 07 | rho | −0.6489 | 1.96E − 09 |
| Income | 0 (0, 0) | 8.56E − 08 | lambda | 0.6118 | 2.38E − 09 |
| 1 Lag | |||||
| Rate ∼ Tobacco + Alcohol + Resin * Herb + LM.Cannabis_x_Resin.THC + Amphetamines + Cocaine + Income | |||||
| Alcohol | 0.08 (0.03, 0.13) | 0.0006 | rho | 0.70286 | 3.41E − 16 |
| Resin | −1.53 (−2.8, −0.26) | 0.0185 | lambda | −0.66616 | 6.38E − 12 |
| LM.Cannabis_x_Resin.THC | 1.8 (0.56, 3.04) | 0.0047 | |||
| Income | 0 (0, 0) | 0.0042 | |||
| 2 Lags | |||||
| Rate ∼ Tobacco + Alcohol + Resin * Herb + LM.Cannabis_x_Resin.THC + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.06 (0.03, 0.09) | 1.09E-05 | rho | −0.73393 | 3.50E − 15 |
| Cocaine | 0.16 (0.02, 0.31) | 0.0296 | lambda | 0.66726 | 1.57E − 13 |
| Income | 0 (0, 0) | 0.0017 | |||
Three geospatial models for hip dysplasia are presented in Table 7. Here terms including daily cannabis use are highly statistically significant from 4.13 × 10−14, 5.39 × 10−8 and 9.04 × 10−15, respectively.
| Parameter values | Model parameters | ||||
|---|---|---|---|---|---|
| Parameter | Estimate (CI) | P-value | Parameter | Value | Significance |
| Additive | |||||
| Rate ∼ Tobacco + Alcohol + Daily.Interpol. + LM.Cannabis + Resin + + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.07 (0.04, 0.1) | 1.58E − 05 | rho | −0.81401 | <2.2E − 16 |
| Alcohol | −0.15 (−0.21, −0.1) | 6.59E − 08 | lambda | 0.544 | 1.92E − 09 |
| Daily.Interpol. | 41.2 (30.52, 51.88) | 4.13E − 14 | |||
| Resin | −1.14 (−2.07, −0.21) | 0.0173 | |||
| Income | 0 (0, 0) | <2.2E − 16 | |||
| Interactive | |||||
| Rate ∼ Tobacco + Alcohol + Daily.Interpol. * Resin + LM.Cannabis_x_Herb.THC_x_Daily.Interpol. + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.05 (0.02, 0.07) | 0.0009 | rho | 0.67333 | <2.2E − 16 |
| Daily.Interpol. | 49.1 (23.03, 75.17) | 0.0002 | lambda | −0.7352 | <2.2E − 16 |
| Resin | −2.07 (−3.7, −0.44) | 0.0128 | |||
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | −11.4 (−15.38, −7.42) | 2.04E − 08 | |||
| Income | 0 (0, 0) | <2.2E − 16 | |||
| Daily.Interpol.: Resin | 356 (227.82, 484.18) | 5.39E − 08 | |||
| 2 Lags | |||||
| Rate ∼ Tobacco + Alcohol + Daily.Interpol. + Resin + LM.Cannabis_x_Herb.THC_x_Daily.Interpol. + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.08 (0.04, 0.11) | 8.63E − 06 | rho | −0.84387 | <2.2E − 16 |
| Alcohol | −0.17 (−0.23, −0.11) | 5.90E − 08 | lambda | 0.55385 | 2.18E − 10 |
| Daily.Interpol. | 47.2 (35.26, 59.14) | 9.04E − 15 | |||
| Income | 0 (0, 0) | <2.2E − 16 | |||
Tables 8 and 9 present final geospatial models for polydactyly and syndactyly, respectively. In all cases, terms including various cannabis metrics remain in final models, demonstrate positive regression coefficients and are highly significant.
| Parameter values | Model parameters | ||||
|---|---|---|---|---|---|
| Parameter | Estimate (CI) | P-value | Parameter | Value | Significance |
| Additive | |||||
| Rate ∼ Tobacco + Resin + LM.Cannabis_x_Resin.THC + Herb + LM.Cannabis_x_Herb.THC_x_Daily.Interpol. + Alcohol + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.05 (0.03, 0.07) | 5.00E − 07 | rho | −0.6055 | 7.84E − 09 |
| Resin | −2.26 (−4.11, −0.41) | 0.017152 | lambda | 0.53977 | 6.06E − 08 |
| LM.Cannabis_x_Resin.THC | 2.72 (1.1, 4.34) | 0.0010 | |||
| Herb | 5.83 (3.56, 8.1) | 4.91E − 07 | |||
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | −4.72 (−6.55, −2.89) | 3.94E − 07 | |||
| Alcohol | 0.08 (0.04, 0.12) | 9.51E − 05 | |||
| Amphetamines | −0.13 (−0.23, −0.03) | 0.0084 | |||
| Cocaine | 0.27 (0.13, 0.41) | 0.0001 | |||
| Income | 0 (0, 0) | 0.0003 | |||
| Interactive | |||||
| Rate ∼ Tobacco + Resin * LM.Cannabis_x_Herb.THC_x_Daily.Interpol. + LM.Cannabis_x_Resin.THC * Herb + Alcohol + Amphetamines + Cocaine + Income | |||||
| Resin | −2.53 (−4.86, −0.21) | 0.0330 | rho | −0.3893 | 0.00817 |
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol., | −6.36 (−4.86, −0.21) | 8.45E − 07 | lambda | 0.3604 | 0.00672 |
| LM.Cannabis_x_Resin.THC | 3.58 (1.51, 5.65) | 0.0007 | |||
| Herb | 5.09 (2.04, 8.15) | 0.0011 | |||
| Alcohol | 0.12 (0.07, 0.18) | 7.29E − 06 | |||
| Amphetamines | −0.17 (−0.29, −0.05) | 0.0064 | |||
| Cocaine | 0.39 (0.23, 0.55) | 1.27E − 06 | |||
| 1 Lag | |||||
| Rate ∼ Tobacco + Resin * LM.Cannabis_x_Herb.THC_x_Daily.Interpol. + LM.Cannabis_x_Resin.THC + Alcohol + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.05 (0.02, 0.07) | 0.0002 | rho | −0.364 | 0.0512 |
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | −8.62 (−12.64, −4.6) | 2.51E-05 | lambda | 0.4529 | 0.00581 |
| Alcohol | 0.09 (0.04, 0.14) | 0.0006 | |||
| Amphetamines | −0.22 (−0.34, −0.1) | 0.0003 | |||
| Cocaine | 0.38 (0.18, 0.58) | 0.0002 | |||
| Income | 0 (0, 0) | 0.0019 | |||
| Resin: LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | 30.8 (14.92, 46.68) | 0.0001 | |||
| 2 Lags | |||||
| Rate ∼ Tobacco * LM.Cannabis_x_Resin.THC + Resin * LM.Cannabis_x_Herb.THC_x_Daily.Interpol. + Herb + Alcohol + Amphetamines + Cocaine + Income | |||||
| Tobacco | 0.06 (0.03, 0.08) | 2.64E − 06 | rho | −0.4738 | 0.00329 |
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | −8.06 (−12.23, −3.89) | 0.0002 | lambda | 0.4033 | 0.00683 |
| Herb | 3.06 (0.43, 5.69) | 0.0226 | |||
| Alcohol | 0.07 (0.02, 0.12) | 0.0038 | |||
| Amphetamines | −0.22 (−0.33, −0.1) | 0.0002 | |||
| Cocaine | 0.34 (0.15, 0.53) | 0.0004 | |||
| Income | 0 (0, 0) | 0.0033 | |||
| Resin: LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | 30.8 (11.34, 50.26) | 0.0019 | |||
| Parameter values | Model parameters | ||||
|---|---|---|---|---|---|
| Parameter | Estimate (CI) | P-value | Parameter | Value | Significance |
| Additive | |||||
| Rate ∼ Tobacco + LM.Cannabis_x_Herb.THC + Herb + Resin + LM.Cannabis_x_Resin.THC + Alcohol + Amphetamines + Cocaine + Income | |||||
| Herb | 3.4 (1.52, 5.28) | 0.0004 | rho | −0.7663 | <2.2E − 16 |
| Income | 0 (0, 0) | 1.45E-05 | lambda | 0.69585 | <2.2E − 16 |
| Interactive | |||||
| Rate ∼ Tobacco * Resin + Herb * LM.Cannabis_x_Resin.THC + LM.Cannabis_x_Herb.THC + Alcohol + Amphetamines + Cocaine + Income | |||||
| Tobacco | −0.07 (−0.1, −0.03) | 0.0004 | rho | −0.77296 | <2.2E − 16 |
| Resin | −8.73 (−12.96, −4.5) | 5.42E − 05 | lambda | 0.6701 | <2.2E − 16 |
| Herb | 3.68 (1.29, 6.07) | 0.0026 | |||
| Income | 0 (0, 0) | 0.0113 | |||
| Tobacco:resin | 0.32 (0.15, 0.48) | 0.0002 | |||
| 1 Lag | |||||
| Rate ∼ Tobacco * Herb + Resin + LM.Cannabis_x_Resin.THC_x_Daily.Interpol. * LM.Cannabis_x_Herb.THC + Alcohol + Amphetamines + Cocaine + Income | |||||
| Income | 0 (0, 0) | 0.0153 | rho | −0.75436 | <2.2E − 16 |
| Tobacco:herb | 0.46 (0.02, 0.9) | 0.0433 | lambda | 0.71334 | <2.2E − 16 |
| 2 Lags | |||||
| LAnomRt ∼ lag(TobRt,2) * lag(pmHerb,2) + lag(Resin,2) + lag(pmResin,2) * lag(Herb, 2) + lag(AlcRt,2) + lag(LAmph,2) + lag(LCocc,2) + lag(MHY,2) | |||||
| Resin | −4.93 (−7.52, −2.34) | 0.0002 | rho | −0.83404 | <2.2E − 16 |
| LM.Cannabis_x_Resin.THC | 3.39 (0.97, 5.8) | 0.0060 | lambda | 0.60161 | 4.43E − 10 |
| Herb | 10.97 (7.63, 14.31) | 1.20E − 10 | |||
| Tobacco: LM.Cannabis_x_Herb.THC | −0.46 (−0.74, −0.18) | 0.0015 | |||
E-values can be extracted from each of these regression models. E-values applicable to panel and geospatial regression models are presented in Tables 10 and 11, respectively. These E-values are listed together consecutively in Table 12. In this table they are separated from their respective regression terms so that their overall magnitude can be compared across all models as a whole corpus of data. It is found that 49 of the 57 E-value estimates are greater than 9 and thus high [102], and all 57 (100%) exceed the threshold for causality of 1.25 [101]. Considering the 95% lower bound of the E-values 40/57 (70.2%) are greater than 9 and 56 (98.2%) exceed 1.25.
| Anomaly | Model and term | P-value | E-value estimate | Lower bound E-value |
|---|---|---|---|---|
| Limbs | Additive | |||
| Herb | 0.00746 | 2049.60 | 13.73 | |
| Interactive | ||||
| Tobacco:herb | 5.14E − 05 | 5.35 | 3.92 | |
| 2 Lags | ||||
| Tobacco:herb | 0.0025 | 4.22 | 2.11 | |
| Limb reduction | Additive | |||
| Resin | 1.12E − 05 | 518.03 | 49.56 | |
| Interactive | ||||
| Tobacco:herb:resin | 9.89E − 05 | 7.06 | 3.41 | |
| 1 Lag | ||||
| Tobacco:herb | 0.0047 | 4.11 | 1.97 | |
| Polynomial in | Quadratic in resin | |||
| Limb reduction | I(Resin^2) | 2.28E − 05 | 2.79E + 07 | 2.29E + 04 |
| Tobacco:herb | 0.0267 | 3.99 | 1.48 | |
| Quadratic—Linear interaction | ||||
| (Resin^2): LM.Cannabis_x_Resin.THC | 1.70E − 05 | 1.34E + 72 | 4.90E + 41 | |
| Quadratic—Quadratic interaction | ||||
| (LM.Cannabis_x_Resin.THC^2): I(Resin^2) | 0.0001 | 3.40E + 65 | 6.14E + 33 | |
| Quartic in resin | ||||
| (Resin^4) | 1.69E-05 | 2.78E + 85 | 1.84E + 49 | |
| Hip | Interactive | |||
| Tobacco:herb | 0.037 | 5.74 | 1.41 | |
| 2 Lags | ||||
| Tobacco:herb | 0.0010 | 19.67 | 7.19 | |
| Polydactyly | Additive | |||
| LM.Cannabis_x_Resin.THC | 6.75E − 08 | 2.37E + 04 | 978.47 | |
| Herb | 0.0001 | 1.11E + 07 | 6.22E + 03 | |
| Interactive | ||||
| Herb | 6.72E − 06 | 1.68E + 28 | 6.39E + 16 | |
| LM.Cannabis_x_Resin.THC | 2.22E − 07 | 2.70E + 05 | 4.39E + 03 | |
| Tobacco: Herb: LM.Cannabis_x_Resin.THC | 6.20E − 05 | 2.02E + 07 | 1.11E + 04 | |
| 1 Lag | ||||
| LM.Cannabis_x_Resin.THC | 6.03E − 05 | 2.72E + 04 | 332.97 | |
| Herb | 0.0030 | 1.21E + 06 | 242.35 | |
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol.: LM.Cannabis_x_Resin.THC | 1.06E-05 | 9.30E + 32 | 2.26E + 19 | |
| 2 Lags | ||||
| LM.Cannabis_x_Resin.THC | 4.60E − 06 | 5.79E + 06 | 1.63E + 04 | |
| Herb | 0.014 | 6.18E + 05 | 33.19 | |
| LM.Cannabis_x_Herb.THC_x_Daily.Interpol.: LM.Cannabis_x_Resin.THC | 0.0083 | 1.98E + 22 | 2.86E + 06 | |
| Resin: LM.Cannabis_x_Herb.THC_x_Daily.Interpol.: LM.Cannabis_x_Resin.THC | 1.65E − 07 | 8.13E + 137 | 3.09E + 91 | |
| Syndactyly | Additive | |||
| Resin | 0.0017 | 87.92 | 8.41 | |
| Herb | 0.0019 | 5.03E + 05 | 247.29 | |
| Interactive | ||||
| Resin | 0.0195 | 72.74 | 3.19 | |
| Herb | 0.0013 | 4.49E + 06 | 823.57 | |
| LM.Cannabis_x_Herb.THC: LM.Cannabis_x_Resin.THC_x_Daily.Interpol. | 0.0256 | 9.15 | 1.79 | |
| 1 Lag additive | ||||
| Resin | 9.58E − 07 | 100.38 | 23.07 | |
| 1 Lag interactive | ||||
| LM.Cannabis_x_Resin.THC | 9.06E − 06 | 4.90E + 31 | 5.28E + 18 | |
| Herb | 0.0003 | 3.93E + 05 | 660.46 | |
| Tobacco: Resin: LM.Cannabis_x_Resin.THC | 2.18E − 05 | 8.34E + 04 | 822.47 | |
| 2 Lags | ||||
| Resin | 0.0147 | 15.44 | 2.51 |
| Anomaly | Parameter | P-value | E-value estimate | Lower bound E-value |
|---|---|---|---|---|
| Limb | Additive | |||
| Herb | 0.0364 | 317.12 | 2.13 | |
| Resin | 0.0006 | 37.33 | 6.57 | |
| Interactive | ||||
| Herb | 0.0047 | 367.64 | 9.43 | |
| 2 Lags | ||||
| Herb | 0.0066 | 837.86 | 10.33 | |
| Limb reduction | 1 Lag | |||
| LM.Cannabis_x_Resin.THC | 0.0047 | 1.14E + 02 | 6.40E + 00 | |
| Hip dysplasia | Additive | |||
| Daily.Interpol. | 4.13E − 14 | 1.09E + 24 | 7.85E + 17 | |
| Interactive | ||||
| Daily.Interpol. | 0.0002 | 1.57E + 31 | 6.17E + 14 | |
| Daily.Interpol.:resin | 5.39E − 08 | Infinity | 4.58E + 143 | |
| 2 Lags | ||||
| Daily.Interpol. | 9.04E − 15 | 6.22E + 26 | 1.28E + 20 | |
| Polydactyly | Additive | |||
| LM.Cannabis_x_Resin.THC | 0.0010 | 1.12E + 03 | 25.41 | |
| Herb | 4.91E − 07 | 1.58E + 06 | 5.06E + 03 | |
| Interactive | ||||
| LM.Cannabis_x_Resin.THC | 0.0007 | 2.52E + 03 | 40.75 | |
| Herb | 0.0011 | 5.19E + 04 | 117.31 | |
| 1 Lag | ||||
| Resin: LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | 0.0001 | 3.35E + 28 | 9.79E + 13 | |
| 2 Lags | ||||
| Herb | 0.0226 | 2.17E + 03 | 4.86 | |
| Resin: LM.Cannabis_x_Herb.THC_x_Daily.Interpol. | 0.0019 | 8.46E + 30 | 4.35E + 11 | |
| Syndactyly | ||||
| Additive | ||||
| Herb | 0.0004 | 4.78E + 03 | 65.42 | |
| Interactive | ||||
| Herb | 0.0026 | 1.43E + 03 | 19.62 | |
| Tobacco:resin | 0.0002 | 2.91 | 1.94 | |
| 1 Lag | ||||
| Tobacco:herb | 0.0433 | 4.99 | 1.21 | |
| 2 Lags | ||||
| LM.Cannabis_x_Resin.THC | 0.0060 | 4.56E + 03 | 17.97 | |
| Herb | 1.20E − 10 | 1.45E + 11 | 7.29E + 07 |
| No. | E-value estimate | Lower bound E-value |
|---|---|---|
| 1 | Infinity | 4.58E + 143 |
| 2 | 8.13E + 137 | 3.09E + 91 |
| 3 | 2.78E + 85 | 1.84E + 49 |
| 4 | 1.34E + 72 | 4.90E + 41 |
| 5 | 3.40E + 65 | 6.14E + 33 |
| 6 | 9.30E + 32 | 1.28E + 20 |
| 7 | 4.90E + 31 | 2.26E + 19 |
| 8 | 1.57E + 31 | 5.28E + 18 |
| 9 | 8.46E + 30 | 7.85E + 17 |
| 10 | 3.35E + 28 | 6.39E + 16 |
| 11 | 1.68E + 28 | 6.17E + 14 |
| 12 | 6.22E + 26 | 9.79E + 13 |
| 13 | 1.09E + 24 | 4.35E + 11 |
| 14 | 1.98E + 22 | 7.29E + 07 |
| 15 | 1.45E + 11 | 2.86E + 06 |
| 16 | 2.79E + 07 | 2.29E + 04 |
| 17 | 2.02E + 07 | 1.63E + 04 |
| 18 | 1.11E + 07 | 1.11E + 04 |
| 19 | 5.79E + 06 | 6.22E + 03 |
| 20 | 4.49E + 06 | 5.06E + 03 |
| 21 | 1.58E + 06 | 4.39E + 03 |
| 22 | 1.21E + 06 | 978.47 |
| 23 | 6.18E + 05 | 823.57 |
| 24 | 5.03E + 05 | 822.47 |
| 25 | 3.93E + 05 | 660.46 |
| 26 | 2.70E + 05 | 332.97 |
| 27 | 8.34E + 04 | 247.29 |
| 28 | 5.19E + 04 | 242.35 |
| 29 | 2.72E + 04 | 117.31 |
| 30 | 2.37E + 04 | 65.42 |
| 31 | 4.78E + 03 | 49.56 |
| 32 | 4.56E + 03 | 40.75 |
| 33 | 2.52E + 03 | 33.19 |
| 34 | 2.17E + 03 | 25.41 |
| 35 | 2049.60 | 23.07 |
| 36 | 1.43E + 03 | 19.62 |
| 37 | 1.12E + 03 | 17.97 |
| 38 | 837.86 | 13.73 |
| 39 | 518.03 | 10.33 |
| 40 | 367.64 | 9.43 |
| 41 | 317.12 | 8.41 |
| 42 | 113.89 | 7.19 |
| 43 | 100.38 | 6.57 |
| 44 | 87.92 | 6.40 |
| 45 | 72.74 | 4.86 |
| 46 | 37.33 | 3.92 |
| 47 | 19.67 | 3.41 |
| 48 | 15.44 | 3.19 |
| 49 | 9.15 | 2.51 |
| 50 | 7.06 | 2.13 |
| 51 | 5.74 | 2.11 |
| 52 | 5.35 | 1.97 |
| 53 | 4.99 | 1.94 |
| 54 | 4.22 | 1.79 |
| 55 | 4.11 | 1.48 |
| 56 | 3.99 | 1.41 |
| 57 | 2.91 | 1.21 |
| Anomaly | No. of terms | Median E-value estimate | Range E-value estimate | Median lower E-value | Range lower E-value |
|---|---|---|---|---|---|
| Limb reduction polynomials | 5 | 3.40E + 65 | 3.99–2.78E + 85 | 6.14E + 33 | 1.48–1.84E + 49 |
| Hip | 6 | 1.09E + 24 | 5.74–Infinity | 6.17E + 14 | 1.41–4.58E + 143 |
| Polydactyly | 19 | 1.58E + 06 | 1.12E + 03–8.13E + 137 | 6.22E + 03 | 4.86–3.09E + 91 |
| Syndactyly | 16 | 4.78E + 03 | 2.91–4.90E + 31 | 65.42 | 1.21–5.28E + 18 |
| Limb | 7 | 317.12 | 4.22–2049.60 | 6.57 | 2.11–13.73 |
| Limb reduction | 4 | 7.06 | 4.11–518.03 | 3.41 | 1.97–49.56 |
Supplementary Table ST19 provides details of the model type, the anomaly, the model structure and the implicated term for the E-values. The table is ordered by anomaly. This table is summarized for both the E-value estimates and their lower bounds in Table 14.
| Term | No. of terms | Median E-value estimate | Range E-value estimate | Median lower E-value | Range lower E-value |
|---|---|---|---|---|---|
| Daily interpolated | 10 | 1.09E + 24 | 9.15–Infinity | 6.17E + 14 | 1.79–4.58E + 143 |
| Resin | 22 | 2.72E + 04 | 2.91–2.78E + 85 | 822.47 | 1.94–1.84E + 49 |
| Herb | 25 | 2049.60 | 4.99–1.68E + 28 | 13.73 | 1.21–6.39E + 16 |
Supplementary Table ST20 plays a similar role to Supplementary Table ST19 but is this time ordered by the regression term. This facilitates a summary and analysis of this table by term. Table 15 presents a summary of this table simplified by the three main cannabis metric terms into daily interpolated cannabis use, resin and herb THC concentration and shows median and ranges for the E-value estimates and lower bounds. The table is ordered in descending order of the E-value estimate. Clearly the daily interpolated cannabis exposure is the most powerful of the three metrics of exposure.
| Comparison | W-statistic | Alternative | P-value |
|---|---|---|---|
| Lower E-value, Daily_v_Herb | 223 | two.sided | 1.24E − 04 |
| Lower E-value, Daily_v_Resin | 172 | two.sided | 0.0106 |
| Lower E-value, Herb_v_Resin | 186 | two.sided | 0.0587 |
| E-value estimate, Daily_v_Herb | 229 | two.sided | 3.36E − 05 |
| E-value estimate, Daily_v_Resin | 172 | two.sided | 0.0106 |
| E-value estimate, Herb_v_Resin | 224 | two.sided | 0.2846 |
The three cannabis metrics are then compared formally using the Wilcoxson test in Table 15. The results shown in this table clearly demonstrate that daily cannabis use interpolated is the most powerful metric of cannabis exposure for both E-value and its lower bound compared to cannabis herb (P = 3.36 × 10−5 and 1.24 × 10−4) and resin (both P = 0.0106) THC concentrations. The comparisons between herb and resin are not significant.
Discussion
Main results
Study data indicate that congenital limb anomalies are strongly positively associated with various metrics of cannabis exposure in Europe in continuous and categorical bivariate analyses and are robust to causal inferential and geospatial multivariable model adjustment for all five groups of limb CAs studied in detail. Data therefore confirm earlier preclinical reports and epidemiological studies linking cannabis and congenital limb anomalies from Hawaii, the USA and earlier and more elementary European analyses [10, 11, 13].
While the present results are based on population cannabis exposure it is important to note that it has been shown repeatedly that maternal cannabis exposure closely tracks rates of community cannabis use [69–76].
Interestingly, a companion paper to this one has shown that the VACTERL syndrome (vertebral, anorectal, cardiac, tracheo-esophageal fistulae/esophageal atresia, renal and limb anomalies) was strongly and causally linked with European cannabinoid exposure (manuscript submitted). Since limb anomalies are part of the VACTERL syndrome this finding also provides evidence from this source confirming the present findings.
Detailed results
Rates of limb anomalies, hip dysplasia and the whole group of limb anomalies were much higher in nations with increasing indices of daily cannabis use (Fig. 3). mEV in the high range exceeding nine on bivariate analysis included hip dysplasia, limb anomalies, talipes and polydactyly (Table 1). In an additive inverse probability weighted panel models the order of teratogenicity was Resin > Herb ≫ Tobacco > Alcohol (Table 2).
Eighty-six percentage of E-value estimates and 70.2% of mEVs of 57 E-value pairs from inverse probability weighted panel models and from spatial models exceeded nine and were thus in the high zone (Table 11). As judged by the mEV the degree of association with metrics of cannabis exposure was hip dysplasia > polydactyly > syndactyly > limb anomalies > limb reductions with median E-value estimates of 3.40 × 1065, 1.09 × 1024, 1.58 × 106, 4.78 × 103, 317.2 and 7.06 and median mEVs 6.14 × 1033, 6.17 × 1014, 6.22 × 103, 65.42, 6.57 and 3.41, respectively (Table 13). When polynomial terms were introduced for cannabis exposure for limb reduction this became the most powerful model with E-value estimate 6.40 × 1065 and mEV 6.14 × 1043. The signal for limb reduction was present at 1 year but disappeared at 2 years in both spatial and panel models. Daily cannabis use interpolated was a more powerful metric of cannabis exposure than herb or resin THC exposure, confirming the important role of intensity of cannabis exposure and presumably related total cannabinoid dose.
The observations from France, Germany, Switzerland and Australia cited in the Introduction point to the suggestion that liberal cannabis policies tend to be linked with cannabinoid food chain contamination which is in turn linked with greatly elevated rates of congenital limb anomalies. Data on cannabis legalization and policies on cannabis in the food chain are not available to this research team at the time of writing. Clearly it would be important to repeat the analyses when these relevant datasets become publicly available.
Choice of anomalies
As the association of limb anomalies with cannabis use had been identified as major issues in preclinical animal models and in both the Hawaiian and USA epidemiological data, investigators were naturally interested to learn how this association would perform in the European dataset where very rich group data were available for analysis. We also chose to study limb anomalies as they are among the most dramatic of CAs and one which stirs considerable public interest. Hip dysplasia was chosen as it obviously had a strong positive relationship with many cannabis metrics on bivariate analysis. Polydactyly and syndactyly were chosen for further study because they had been identified in the Hawaiian series and seemed to allow us further opportunity to interrogate the relationship between prenatal cannabis exposure and disordered limb development in addition to demonstrating an obviously strong positive relationship with cannabis resin THC concentration in Fig. 2.
Qualitative causal inference
In 1965, Hill elaborated nine criteria as the bases for causal inference in epidemiology. It is noteworthy that the present findings fulfill all of these criteria: strength of association, consistency among studies, specificity, temporality, coherence with known data, biological plausibility (note mechanistic comments below), dose–response curve, analogy with similar situations elsewhere (especially the USA and Hawaii) and experimental confirmation [105].
Quantitative causal inference
The use of inverse probability weighting in all panel models has been shown to transfer the analysis from the purely observational context to a pseudorandomized context from which it is appropriate to draw causal inferences [106]. The use of E-values quantifies the degree of association required of an extraneous covariate which has not been controlled in the present analysis with both the outcome of interest and the exposure of concern to void an apparently causal association. The finding therefore in the present work that 86.0% of E-value estimates and 70.2% of mEVs were in the high range clearly places the findings in a position which are robust to extraneous or uncontrolled confounding and provides the reader with a high degree of confidence in the reported results. Application of these two key techniques of formal quantitative inference in the present study has the effect of strongly reinforcing study outcomes in a robustly causal inferential framework.
Mechanisms
Morphogen gradients
Limb development commences from Days 22 to 24 of human gestation for the arm and from Days 24 to 26 for the leg. Like many tissues limb morphogenesis is controlled by 3D mutually antagonistic morphogen gradients in each of the three axes proximal to distal, anterior to posterior and left and right lateral [14]. Gradients between retinoic acid proximally antagonize gradients of fibroblast growth factor 4 (FGF4) and Wnt coming from the distal advancing tip. The limb grows in length under the influence of FGF released from the advancing anterior epidermal ridge. The upper, middle and lower parts of each limb form under the influence of Meis homeobox 1/2 (Meis1/2), Homeobox A11 (Hoxa11) and Hoxa13 respectively. The limb stops growing when the distance from the sonic hedgehog source proximally becomes excessive. Gradients of sonic hedgehog (shh) posteriorly antagonize gremlin gradients from anterior tissues. Digital development and finger growth happens under the control of FGF8 and apoptotic cell death in the intervening interdigital rays is controlled by retinoic acid receptor B, and bone morphogenetic proteins (BMPs) BMP2, BMP4 and BMP7 and Msh homeobox 1 (Msx1). A sonic hedgehog source medially is transduced by expression of its receptor patched in the medial tissues and the signal dissipates more laterally [107]. Wave-like intersecting gradients of shh and Msx1 control the growth of digits and the loss of the intervening tissue in the web spaces.
Cannabinoid inhibition of morphogens
Many of these morphogen gradients are disrupted by cannabinoids including BMPs [108–110], retinoic acid [111–113], Wnt signaling [114–119], FGF [120, 121] and sonic hedgehog [122]. It therefore becomes easy to see how gradient disruption at critical periods might interfere with this delicately balanced and finely tuned process.
Epigenomic controls
One recent carefully designed study identified 163 differentially methylated regions (DMRs) of DNA in human and rat sperm during cannabis dependence and 127 DMRs 11 weeks after confirmed cannabis abstinence [45]. Eleven weeks is the period of one sperm cycle in the human male. These DMRs affected hundreds of genes.
Five genes were identified during the period of cannabis abstinence which together affected functions such as limb morphogenesis and development, skeletal development and apoptosis of limb bud cells (page 354, P = 0.005470). These five genes were BMP4, CHD7 (Chromodomain helicase DNA binding protein 7), GLI3 (Gli zinc family finger 3), MEGF8 (Multiple EGF-like domains 8) and TMEM107 (Transmembrane protein 107) [45].
As noted above BMP4 is the natural antagonist of sonic hedgehog [14]. GLI3 is involved in transducing signals through part of the sonic hedgehog pathway [123, 124]. MEGF8 and TMEM107 are also involved with sonic hedgehog signaling [125–127].
GLI mutations are associated with anomalies of the head and brain and polydactyly and syndactyly (Greig syndrome of cephalopolysyndactyly) [123].
MEGF8 mutations are associated with skull, finger and toe abnormalities including syndactyly, brachydactyly and polydactyly [128]. Six different mutations of MEGF8 have been shown to cause Carpenter syndrome which also includes mental retardation. Carpenter syndrome can also include left right switching of major organs such as the heart (dextrocardia), the great vessels (transposition) or all the organs (situs inversus) [128].
TMEM107 is found in the transitional zone at the base of primary cilia. Both cilia and TMEM107 transduce morphogenic sonic hedgehog signals [127]. Cilia are involved in embryonic patterning left–right specification, neural development, skeletal formation and organogenesis [127]. Ciliopathies include anomalies of fingers, neural patterning defects, microphthalmia and skeletal abnormalities.
BMP4 is critically involved in heart development, microphthalmia, otic placode development and orofacial cleft formation [129]. It is widely expressed in the bladder and controls development of the ureteric bud and ureteric development [130]. It is also involved in limb bud development and development of lungs, face, liver and teeth. It is also involved critically in brain development and neurogenesis from both the subdentate and subventricular zones [131–133]. Lens anomalies and cataract and orofacial clefts also occur in BMP4 mutant animals [134, 135].
Importantly cannabis withdrawal was marked by the significant finding of increased apoptosis of limb bud cells (page 356, one gene, BMP4, P = 0.00701). Clearly this would fit nicely with clinical observations of amelia and phocomelia related to apoptosis of limb bud cells at critical periods of embryogenesis.
Polydactyly was identified in two functional annotations both in cannabis withdrawal as preaxial polydactyly (three genes, GLI3, MEGF8 and TMEM107, P = 0.0000919, page 333) and polydactyly (four genes, BMP4, GLI3, MEGF8 and TMEM107, P = 0.0033 page 350).
Syndactyly was not identified alone but was identified as polysyndactyly in cannabis withdrawal (two genes, GLI3 and MEGF8, P = 0.0017, page 345).
It is clear from this brief review that both direct and epigenetically mediated cannabinoid disruption of these key morphogen gradients during embryonic life can potentially disturb not only limb growth, patterning and development but can also give rise to many other anomalies which have been described in association with prenatal cannabis exposure also.
Exponential genotoxic effects
A rich and varied literature indicates that cannabinoid genotoxicity exhibits an exponential dose–response effect on mitochondrial inhibition, actin, tubulin, RNA, DNA and protein synthesis, chromosomal breaks, cell growth and viability, mutagenesis, epigenomic transmission of addictive behaviors and micronucleus formation which are all related to mutagenic processes [31, 42, 58–60, 62, 122, 136–142].
Therefore, a real concern must exist that the rising triple convergence of cannabis prevalence of use, intensity of daily use and cannabinoid potency in marketed product [49–51, 143] will collide with this strongly exponential genotoxic dose–response curve giving rise to mini-epidemic of CAs such as increased limblessness in northeastern France identified in the current analysis [64–66] or the accelerating rates of atrial septal defects in Kentucky and Mississippi in the USA [144], creating an apparent step-like hike in patterns of clinical teratogenesis.
Strengths and limitations
This study has a number of strengths and limitations. Its strengths include using arguably the largest most comprehensive database of CAs globally, using two forms of multivariate regression, the use of the key tools of causal inference inverse probability weighting and E-values to work with a pseudorandomized dataset and draw causal conclusions, the use of geospatiotemporal modeling, and the use of random forest regression for variable selection. Limitations of this analysis include those of many epidemiological studies in that individual cannabis use and exposure data were not available. Linear interpolation was used to complete the daily cannabis use data as described and needs to be born in mind when considering the reported results. Substance and limb data are whole population estimates and hence cannot illuminate differences in behavioral and genetic sub-populations.
Generalizability
Study findings are derived from arguably the largest most comprehensive dataset globally and analyzed within a pseudo-randomized analytical paradigm that fulfills criteria of epidemiological causality. For these reasons we are confident that study results are generalizable to other contexts, with the proviso that data are of sufficient size and quality to contain the required information. Importantly we also note the consistency of the present results with published study results reported elsewhere from Hawaii and the USA [10, 11].
Moreover, the strong background of biologically plausible explanatory frameworks which explains findings mechanistically strengthens confidence in both findings and generalizability.
Conclusion
In summary, data confirm earlier reports from animals and epidemiological studies in Hawaii and the USA generally linking various metrics of cannabis exposure with a spectrum of limb anomalies using a pseudorandomized causal inferential paradigm and in their native space–time setting. By using the techniques of inverse probability weighting and E-values it has been demonstrated that the present results fulfill epidemiological criteria for causal relationships. Results begin to complete a detailed portrait of cannabis genotoxicity more broadly which also includes aspects of cancer induction and cellular and organismal aging. Great concern is also expressed at the well-demonstrated exponential dose–response relationships of cannabinoid and genotoxicity first studied in the laboratory and now reported in multiple independent epidemiological studies associated with rising levels of cannabinoid exposure in many countries [49–51, 143]. In many of the present analyses cannabis is shown to be a much more potent teratogen than tobacco or alcohol indicating, particularly in light of the known exponential dose–response curve, that access to cannabinoids should be carefully restricted to protect the food chain, to preserve community genomic and epigenomic health, to forestall a series of teratological–mutagenic epidemics, to prevent enhanced carcinogenesis, and to preserve the epigenomic inheritance of coming generations.
Supplementary Material
Acknowledgements
All authors had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.
Supplementary data
Supplementary data are available at EnvEpig online.
Data availability
All data generated or analyzed during this study are included in this published article and its supplementary information files. Data along with the relevant R code have been made publicly available on the Mendeley Database Repository and can be accessed from these URLs: 10.17632/vsmmmkncsd.1 and 10.17632/vd6mt5r5jm.1.
Conflict of interest statement
The authors declare that they have no competing interests.
Funding
No funding was provided for this study. No funding organization played any role in the design and conduct of the study; collection, management, analysis and interpretation of the data; preparation, review or approval of the manuscript; and decision to submit the manuscript for publication.
Ethics approval and consent to participate
The Human Research Ethics Committee of the University of Western Australia provided ethical approval for the study to be undertaken 24 September 2021 (No. RA/4/20/4724).
Consent for publication
Not applicable.