Effects of cannabis on congenital limb anomalies in 14 European nations: A geospatiotemporal and causal inferential study.

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 1:

Paneled bivariate scatterplot of limb CARs against various substances

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.

Figure 2:

Paneled bivariate scatterplot of limb CARs against various metrics of cannabis exposure

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 3:

Paneled bivariate scatterplot of limb CARs over time by grouping of nations into rising or static rates of daily cannabis use. See the Methods section for details of national grouping

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.

Figure 4:

Time sequential map graph of log rates of limb anomalies across Europe for studied nations

Figure 5:

Time sequential map graph of log rates of limb reduction anomalies across Europe for studied nations

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.

Figure 6:

Time sequential map graph of daily cannabis use interpolated across Europe for studied nations

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 7:

Bivariate colorplane sequential map graph series of log rates of limb anomalies by cannabis resin THC concentration across Europe for studied nations. Relative rates are as described in the colorplane key

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.

Figure 8:

Bivariate colorplane sequential map graph series of log rates of limb reduction anomalies by cannabis resin THC concentration across Europe for studied nations. Relative rates are as described in the colorplane key

Figure 9:

Bivariate colorplane sequential map graph series of log rates of polydactyly anomalies by cannabis resin THC concentration across Europe for studied nations. Relative rates are as described in the colorplane key

Figure 10:

Bivariate colorplane sequential map graph series of log rates of syndactyly anomalies by cannabis resin THC concentration across Europe for studied nations. Relative rates are as described in the colorplane key

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).

Figure 11:

Bivariate colorplane sequential map graph series of log rates of hip dysplasia anomalies by cannabis resin THC concentration:last month cannabis use:daily cannabis use interpolated across Europe for studied nations. Relative rates are as described in the colorplane key

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].

Table 1:

Positive Significant Bivariate Slopes from Linear Models for Selected Limb Anomalies

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

Note that the values in this table are ordered in descending order by the lower bound of the E-value.

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.

Figure 12:

Paneled boxplots of various CARs by extreme daily cannabis use quintiles. Nonoverlapping notches indicate statistical significance

Figure 13:

Paneled boxplots of various CARs by extreme quintiles of daily cannabis use interpolated:cannabis resin THC:last month cannabis use. Nonoverlapping notches indicate statistical significance

Figure 14:

Paneled boxplots of various CARs by extreme quintiles of daily cannabis use interpolated:cannabis herb THC:last month cannabis use. Nonoverlapping notches indicate statistical significance

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.

Table 2:

Cannabis-Related Bivariate Categorical Comparisons of Highest and Lowest Monthly Cannabis Use Quintiles

		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

Note that the values are ordered alphabetically by limb anomaly. The anomaly “limb” relates to the whole class of limb anomalies.

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.

Table 3:

Positive Significant Coefficients from Serialized Additive Multivariable Mixed Effects Models in Purrr by Anomaly

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.

Table 4:

Summary of Serialized Additive Mixed Effects Models in Purrr Shown in Table 3 by Substance

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

Figure 15:

Summary data for additive mixed effects models by substance. (A) Number of anomalies, (B) sum of P-value exponents, (C) mean value of P-value exponents, (D) total percentage increase in marginal effect from the minimum to the maximum value of covariate, (E) mean percentage increase in the marginal effect of the increase of the covariate from the minimum to the maximum and (F) median percentage increase in the marginal effect of the increase of the covariate from the minimum to the maximum

Figure 16:

Summary of E-values for additive mixed effects models by substance. (A) Sum of E-value estimates, (B) mean E-value estimates, (C) median E-value estimates, (D) sum of lower E-values, (E) mean of lower E-values and (F) median of lower E-values

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.

Table 5:

Multivariable Geospatial Regression Models for Limb Anomalies

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.

Table 6:

Multivariable Geospatial Regression Models for Limb Reduction Anomalies

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.

Table 7:

Multivariable Geospatial Regression Models for Hip Dysplasia

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.

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.

Table 12:

E-Value Lists for All Values

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

Cannabis-related terms have been chosen for inclusion in this table.

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.

Table 14:

E-Value Summary by Cannabis Metric

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.

Table 15:

Wilcoxon Test Results for Inter-group Comparisons for Main Cannabis Metrics

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.

Table 11:

E-Values from Spatial Models

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

Table 13:

E-Value Summary by Anomaly

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

Cannabis-related terms have been chosen for inclusion in this table.

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.

Click here for additional data file.

Table 8:

Multivariable Geospatial Regression Models for Polydactyly Anomalies

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

Table 9:

Multivariable Geospatial Regression Models for Syndactyly Anomalies

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

Table 10:

E-Values from Panel Models

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