Logit-Normal Mixed Model for Indian Monsoon Precipitation

Describing the nature and variability of Indian monsoon precipitation is a topic of much debate in the current literature. We suggest the use of a generalized linear mixed model (GLMM), specifically, the logit-normal mixed model, to describe the underlying structure of this complex 5 climatic event. Several GLMM algorithms are described and simulations are performed to vet these algorithms before applying them to the Indian precipitation data. The logit-normal model was applied to light, moderate, and extreme rainfall. This work provides a novel use of GLMM and promotes its’ 10 addition to the gamut of tools for analysis in studying climate phenomena.


Introduction
Explanation of Indian monsoon precipitation has been a challenging problem in physics as well as data analysis.In this paper, we focus on statistical analysis of the summer monsoon precipitation data, to provide insight symbiotic with deterministic physics modeling.Previous statistical analysis studies regarding precipitation in Indian monsoons have explored two main areas-identifying methodology of data analysis and covariate selection.
The establishment of appropriate statistical methodology for explanation and prediction of precipitation, while simultaneously capturing underlying variability, is paramount.These methods are used in identification of trends for prediction, however, trends tend to be inconsistent across studies and may relate to linked variability on different temporal and spatial scales as noted by Turner and Annamalai (2012).
For instance, Goswami et al. (2006) used daily central Indian rainfall and found rising trends in frequency and magnitude of extreme rain events along with decreasing light and moderate rainfall.While validating their 2006 study, Ghosh et al. (2012) indicated increasing spatial variability in observed Indian rainfall extremes.They also found that moderate rainfall increased in central India despite a decreasing trend in occurrence of moderate rainfall.For high and extremely high rainfall, they noted a few locations experienced a significant upward or downward trend, however, most grid boxes showed a lack of trend.

40
A similar study conducted by Ghosh et al. (2009) used a finer spatial scale and indicated a mixture of increases and decreases of extreme rainfall events dependent on location.An increasing trend in exceedances of 99th (extreme) percentile daily rainfall was discovered by Krishnamurty et al. 45 (2009).On the other hand, they stated many parts of India exhibited a decreasing trend for exceedances of the 90th (moderate to extreme) percentile.Increases in the frequency of both light and moderate to extreme rainfall events were observed in Singh et al. (2014), along with decreasing proba-50 bility of regional rainfall events and higher variability in the intensity of these events.
These studies utilized parametric -regression, extreme value theory, time series methods -and nonparametric statistical techniques, yet their lack of unanimity suggests impor-55 tant properties of the Indian monsoon remain partially misunderstood.
In view of the above, we propose adding the generalized linear mixed model (GLMM) as a potential framework for analysis of Indian monsoon precipitation data.A GLMM is 60 a broader framework compared to the standard (linear, loglinear, logistic, or other) regression in that there are random effects involved.This implies part of the signal is random, and changes from one set of circumstances to another.In the current context, a GLMM may be suitable for capturing lo-65 cal, instantaneous variability.Such local variability may arise from cloud and other physical micro-properties.When there is no such local variability, an appropriate variance component in the GLMM would be zero, thus, recovering the true underlying "fixed-effects" regression model.

70
The second principal focus of literature has been identifying relevant covariates for study of Indian monsoon precipitation.Certain oscillations are commonly useful predictors for precipitation.For instance, the synoptic activity in- Li and Yanai (1996); Prell and Kutzback (1992); Turner and Annamalai (2012)) is cited as a driver of the monsoon.
Several other climactic predictors of monsoons have been proposed in the literature including Himalayan/Eurasian snow extent (Kumar et al. (1999)), Pacific trade winds (Li and Yanai (1996)), atmospheric CO 2 concentration (Prell and Kutzback (1992)), and tropospheric temperature difference (Xavier et al. (2007)).Unfortunately, none have been conclusively attributed for the monsoon rainfall which suggests an intricate relationship between some or all of these factors.
Because explicit attribution to covariates may not be possible, GLMM is a logical model for Indian monsoon precipitation.It allows underlying randomness to drive observed data in a particular hierarchy while still accounting for hypothesized drivers of rainfall.
This paper provides an introduction on extending GLMMs to climate applications.Three paradigms of estimation-approximate likelihood, method of moments, and Bayesianwere tested using four separate algorithmic implementations.The methods of estimating GLMMs were penalized quasilikelihood, penalized iteratively reweighted least squares, method of simulated moments, and data cloning.The theory and limitations of these estimates are described in detail, then utilized in simulations to test the validity of the methodology.
Simulation findings showed that penalized quasilikelihood was not accurate for the given application, thus, the three remaining methods were used to fit logit-normal models with random intercepts by weather station for Indian summer rainfall data in light, moderate, and extreme rainfall classifications.Maximum temperature and elevation were consistently significant in the models aligning with the physics of precipitation.∆T T -tropospheric temperature difference-was also significant for many of the models.The most meaningful finding was a random effect by weather station was non-negligible in many of the models.This provides further credibility to the methodology in applications to climate.Overall, we feel GLMMs could be a significant addition to data analytics in climate applications.The rest of the paper is organized as follows.Section 2 gives a short background on GLMMs and in particular, elucidates the logit-normal model.The theory of the chosen estimation methods for GLMMs are discussed in Section 3. Section 4 furnishes the results of several simulations using these existing methods.Section 5 applies these methods to monsoon precipitation data from India.Finally, Section 6 presents conclusions and future work in this area.
2 Overview of Generalized Linear Mixed Models

Exponential Families
Before discussing GLMMs, we provide preliminaries on the key component-use of an exponential family for the ob-130 served data.An exponential family probability mass/density function (pmf/pdf) has several unique properties conducive to modeling.For further discussion of these properties, refer to Ch. 2 of Keener (2010).For simplicity, consider a univariate random variable Y distributed as an exponential family.

135
The canonical form of the pmf/pdf then can be written: Note a(•) is a function of a dispersion parameter φ, r(•, •) is a function of data and the dispersion parameter, and c(•) is a 140 function of parameters and is known as the cumulant function.The statistic y is complete and sufficient; it is known as the canonical statistic with corresponding canonical parameter η.An exponential family can be written in a more general fashion compared to (1), but will not be discussed here for 145 simplicity.

Model Description
A GLMM is a probability model with a hierarchical structure.Given the latent unobservable second layer, known as random effects, the top layer has a pdf/pmf following an ex-150 ponential family distribution.Assume the observed data are independent conditional on the random effects, and that we have i = 1, . . ., N observations.Define the i th response as Y i .
Then, we can write a GLMM as: Level 1: Level 2: 160 The p-components of the vector β are called fixed effects.The random effects covariance Σ is a function of the q-dimensional (σ 1 , . . ., σ q ) known as variance components.
Fixed covariates are represented by the p × 1 vector x i .Random covariate vectors for the i th data point and r th vari-165 ance component can be denoted by a m r × 1 vector z ir .We combine the vectors for each variance component to form the random covariate vector z i = z T i1 , . . ., z T iq T of length M = q r=1 m r .The random effects vector, U , follows an M -dimensional normal distribution with mean vector 0 and covariance matrix Σ.
In (2), b(•) is a function of only the canonical parameter η i .The "linear" part of GLMM comes from the fact that η i can be represented as a linear function of the fixed and random parameters.r(•, •) is a function of the data and φ.As before, a(•) is a dispersion function.
To illustrate this form, we consider a commonly used version of this model, the logit-normal GLMM.The randomintercept form of this model is: Level 2: U i ind.
Notice that in this model, z i is a vector with a 1 in the i th position and 0's in all other positions.
Returning to the generic form of the GLMM, the assumption of conditional independence among observations implies the density of Y |u is and that the joint density of (Y, U) is However, since random effects are unobserved, in order to utilize the observed data likelihood, one must find the marginal distribution with respect to the observed data Y only.The log-likelihood is then: This integral is rarely analytically tractable.Thus, maxi-200 mum likelihood estimation (which is usually preferable when possible) is very difficult.Many methods for inference have been proposed.Variants of the most popular methods are examined in §3.
3 Methods for Estimating in GLMM 205

Likelihood Approximation Methods
Both methods discussed in the following sections make use of a technique known as Laplace approximation to approximate an integral by a normal distribution (Tierney and Kadane (1986)).Then (10) can be written as: 215 Then, we can express the log likelihood as follows: (β, Σ; y) = log e h(u) du.
This expression can now be approximated.To use the approximation, one first needs the maximizer of the integrand.

220
Let u 0 be the maximizer of e h (u) .Then a Taylor expansion around u 0 yields the approximation to the log-likelihood, Penalized quasi-likelihood (PQL) was proposed by Breslow and Clayton (1993) to approximate the high-dimensional integral using Laplace approximation as a method for obtaining u 0 and ∂ 2 h(u)/∂u∂u T .Filling in the details of ( 13): This equation is differentiated with respect to u and β respectively.Further approximations are made within the derivatives because Σ is also unknown.The approximate derivatives are used to form estimating equations for the 235 mean parameters.For more detailed discussion of these approximations, please refer to McCulloch and Searle (2010).The same estimating equations arise from jointly maximizing 240 with respect to u and β.These equations are solved by using Fisher scoring as an iterated reweighted least squares (IRLS) problem.The quasilikelihood, log f (y|u, β), is optimized taking into account 245 the penalty, 1 2 u T Σu.This penalty term has a shrinkage effect, i.e. forces values of u to be closer to zero.
Variance components in Σ are subsequently estimated using a restricted maximum likelihood approach.Further details on the estimation algorithm are found in §2 of Breslow 250 and Clayton (1993).
The function in R which computes PQL estimates is glmmPQL{MASS}.PQL is reasonably accurate when data are approximately normal and can be very fast.However, Lin and Breslow (1996) and others have criticized this method 255 for it's bias in highly non-normal data.It is especially bad in binomial data with a small sample size or true probabilities near zero or one.Reliance on the quadratic expansion of the log-likelihood is appropriate with normal random effects, yet it is very difficult to assess normality of these unobserved 260 effects.

Penalized Iteratively Reweighted Least Squares
Another approach to likelihood approximation is presented by Bates (2010).The main difference from PQL is that it attempts to approximate the true likelihood rather than the quasi-likelihood.
To understand the approach, first, let U ∼ N M (0, Σ).Consider the decomposition of the random effects covariance matrix Σ = ΓΓ T .Then, U = ΓV where V ∼ N M (0, I m ).This implies that the canonical parameter in (3) can be written as Substituting in v to (9), we note that f (y, v) is proportional to f (v|y).Thus, v 0 is found to maximize f (v|y).
The penalized iteratively reweighted least squares (PIRLS) algorithm is as follows.
1. Given starting values for β, Σ, and v 0 , evaluate η, µ Y |v , and 3. Update the weights, W , and check for convergence.If not converged, go to step 2.
Once the conditional mode ṽ is determined, a Laplace approximation to the deviance (-2* log-likelihood) is evaluated at ṽ.This evaluation may alternatively be done by Gauss-Hermite quadrature which is discussed further in Bolker et al. (2009).The function in R used to compute estimates is glmer{lme4} .This method can experience similar problems to PQL in cases where the random effects are nonnormal.Gauss-Hermite quadrature can allay some of these issues, but is only computationally feasible for small numbers of random effects.

Method of Simulated Moments
Jiang (1998) describes methodology known as the method of simulated moments (MSIM).The method first derives a set of sufficient statistics.Estimating equations are then obtained by equating sample moments of sufficient statistics to their expectations.
Referring to the model elucidated in ( 5)-( 7), let the dispersion function be a(φ) = wi φ where w i is a weight depending on the exponential family of the response.Let θ = (β, σ 1 , . . ., σ q ).Restrict all elements of the z i to be either 0 or 1. Represent where V r ∼ N mr (0, I mr ).Then, where C(•, •, •) represents the other portion of the function.This yields canonical parameters (β/φ, σ 1 /φ, . . ., σ q /φ) 315 with corresponding sufficient statistics 320 Note that the expectations on the right hand side are functions of the parameters while the formulae on the left hand sides are functions of data only.Since the expectations are not available, they must be estimated by Monte Carlo sim-325 ulation.The system of equations can then be solved for the parameters using the Newton-Raphson algorithm.
We implemented this method in a newly created R program.As shown in Jiang (1998), this method is consistent and is potentially computationally less intensive than a 330 Markov Chain Monte Carlo (MCMC) method.

Data Cloning
GLMM estimates can be produced in a traditional Bayesian framework; one must choose priors for the parameters of interest and calculate the posterior distribution by multiplying 335 the prior densities by the likelihood, L(β, Σ|Y ), corresponding (10).One may then use MCMC to generate a dependent sample from the posterior distribution from which estimates can be derived based on strong laws.Lele et al. (2010) derived a method called data cloning to 340 be used in conjunction with MCMC.The algorithm can be summarized in the following three steps.First, create a kcloned data set y k =(y,y,...,y) where the observed data vector is repeated k times.Choose a prior distribution π(β, Σ).
Then, the posterior distribution, π k (β, Σ|Y ), which corre- Under regularity conditions as k → ∞, where (β, Σ) is the MLE of (β, Σ) and S is the Fisher information matrix of the original data.Thus, large k means the posterior distribution is nearly degenerate at the MLE.
To generate a dependent sample from the posterior distribution π k (β, Σ|Y ), one may use an appropriate MCMC algorithm, such as a Gibbs sampler or Metropolis-Hastings algorithm.
Finally, one can calculate the sample means and variances of the components of (β, Σ).Estimates of MLEs for (β, Σ) 360 correspond to these sample means and approximate variances of estimated MLEs correspond to k times the posterior variance of the original data as seen in ( 25).
This method was implemented using dclone{dclone} discussed in Solymos (2010) which relies on the well-reputed 365 BUGS language for estimation of hierarchical models.The method is computationally intensive, and it may prove difficult to assess convergence as with any MCMC implementation.
To quantitatively describe the estimation discrepancy between µ and μm,n , we used squared error loss, Because squared error loss is criticized for a bounded parameter space, we used Stein's loss, to measure how well σ 2 was estimated by σ2 m,n .A combined loss was then calculated as Ideally, as m, n → ∞, G(μ m,n , σ2 m,n ) → 0.

Simulation Estimation Analysis
The estimation results are displayed in Tables 1-4.All methods failed to reasonably estimate both µ and σ 2 in the smallest scenario with 10 subjects and 2 observations each.

400
This was expected because there are not enough replications within subject to get a meaningful estimate of a variance by subject.
All other settings for MSIM, dclone, and glmer estimated µ within 2 standard errors.These methods also pro-405 vided reasonable estimates of σ 2 for settings other than those with 10 subjects.The combination of the loss for the two estimates went to 0 quickly for all three methods.In general, estimation by these three methods were unbiased.
The method glmmPQL did not converge to the true values 410 of (µ, σ 2 ) as evidenced by combined loss greater than 0 for all settings.Further, this method displayed an underestimating bias in both parameters.Also, this function in R could not produce estimates for some of the 100 data sets in each setting. 415

Simulation Speeds
Subsequently, we tested 4 of the 16 simulation settings to determine computing speed of the estimation methods.The settings used were combinations of (50, 200) subjects with (10, 200) observations.The system.time() command in 420 R was used to record times.Simulations were independently run on 4 computers, and each estimation method was tested in sequence in one R script on a single core.Computer specifications can be found in the appendix.We implemented MSIM in two ways for the speed test.

425
In the intercept-only model ( 26)-( 28), it is possible to use a simple algorithm for estimation.However, a more general form of the algorithm is needed for problems including fixed covariates.This form relies on matrices and does not work with large data sets at this time.These methods are referred 430 to MSIM Fast and MSIM Slow, respectively.Results were similar for each of the 4 computers, therefore, only one of the sets of results are shown in Table 5.The results indicated that glmmPQL was fastest in the 50 subject cases and glmer was fastest in the 200 subject cases.These 435 two likelihood methods were the fastest due to the nature of the approximations that they make.The Bayesian method, dclone, was slower at about 4 to 25 minutes to produce estimates.The simple algorithm of MSIM Fast was faster than dclone and slightly slower than approximation meth-440 ods taking 3 to 6 seconds per run.The MSIM Slow method was much slower ranging from 1.5 minutes to nearly 4 hours.The case with 200 subjects 200 observations could not be handled by this method because the size of matrices and vectors exceeded the storage capacity allowed by R.

6
L. Dietz and S. Chatterjee: Logit-Normal Mixed Model for Indian Monsoon Rainfall Precipitation

Simulation Conclusions
In general, glmer and glmmPQL were the fastest, but glmmPQL did not provide good estimates.Parameters were accurately estimated by dclone, but this was much slower than the approximation methods.
In the intercept-only implementation, MSIM proved to be fast, however, was much slower than other methods in it's matrix version and failed when too many observations were used.It should be noted that tuning parameters within each of the methods, such as convergence criteria for MSIM or number of MCMC samples in dclone may impact computing time significantly.
Based on the output of these simulations, the glmmPQL method was not consistent.The other three methodsglmer, dclone, and MSIMprovided estimates with reasonable accuracy.Therefore, these estimation methods are used to fit models for Indian monsoon precipitation data.

Data
Multiple datasets have been used to study Indian monsoon precipitation in the literature of different temporal and spatial granularity.However, the initial goal was to develop and test the methods on widely and freely available datasets for the purpose of understanding the usefulness of GLMM in this context.This led to our selection of the data sources described below.
We chose the National Climatic Data Center (NCDC) 1 in the National Oceanic and Atmospheric Administration (NOAA) to gather latitude ( • ), longitude ( • ), elevation (meters), and daily minimum and maximum temperatures ( • C).These data were collected from 1 January 1973 to 31 December 2013.Data were queried for all available Indian stations in the database.This data source was developed for a wide variety of potential applications, including climate analysis and monitoring studies that require data at a daily time resolution.Quality assurance checks are routinely applied to the full dataset according to Menne et al. (2012).
We note this data had a large amount of missing observations, therefore, only stations with at least 5 observations were included in analysis.One year in particular, 1975, did not contain enough data to be included in the analysis.To elucidate this missingness, on 25 August 2012, there were 33 stations with missing (NA) values, 12 stations with precipitation of 0 mm, and 31 stations with greater than 0mm precipitation.This implies several stations were not included in the data for this day and in general, stations included change over time.Fig. 1 illuminates the rainfall on this date.
We also included several other covariates of interest.The first was tropospheric temperature difference (∆T T ): the air 1 http://www.ncdc.noaa.gov/temperature averaged between the levels 600 hPa and 200 495 hPa.The hypothesis that Indian ocean warming leads to reduction in ∆T T which in turn reduces monsoon circulation is noted in Xavier et al. (2007).Thus, the inclusion of this covariate in the models was relevant.Data were collected from the National Centers for Environmental Prediction (NCEP) 500 Reanalysis site2 .
As stated in Wang ( 2006) and other literature, Indian rainfall is strongly associated with ENSO, and onset of discharge in Nino-3.4region can lead to drought in India.The occurrences of precipitation extremes are thought to be fewer in 505 drought years.The Nino-3.4 monthly anomaly series was gathered for inclusion in the models from the NCEP site sponsored by NOAA3 .
The Indian Ocean Dipole (IOD) is an irregular oscillation occurring in the Indian Ocean.It is commonly measured by 510 the Indian Dipole Mode Index (DMI) which takes the difference between SST anomalies in the western and eastern Indian Ocean.Non-ENSO drought years are associated with DMI thus, this is a relevant covariate for inclusion in modeling.This index was only available for 1973-2010 models and 515 data were procured from the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) site 4 .
We note that analysis of monsoon precipitation using thresholds was previously done in Krishnamurty et al. (2009).Rather than use a fixed threshold for the entirety of 520 India, they utilized data derived percentile thresholds which changed depending on spatial location.However, their research was focused on trend analysis.Since we were able to include spatial covariates, we only consider fixed thresholds for the entire country found in Attri and Tyagi (2010).
All stations were were marked each day with indicators of these categories to be used in the modeling.Only obser-530 vations considered to be within monsoon season were used.This conservatively included the time period from May 1 to October 31 (184 days) for each year.We fit models for each year (excluding 1975) from 1973-2013.To account for spatial variability, we fit a random intercept by weather station 535 (WS) in the following logit-normal model: Level 2: U W S ind.

Results of GLMMs
To aid interpretation and provide a basis for comparison among models, we performed tests of significance for both fixed and random parameter estimates.We also give results from a goodness-of-fit test for each year's model.
In order to provide tests of significance for fixed effect coefficients, we propose the following procedure: 1. Run a generalized linear model (GLM) with all eligible fixed covariates.
2. Run a GLM with all eligible fixed covariates except the one we are testing.
3. Do a likelihood ratio test (LRT) to compare these models and get a p-value from the asymptotic χ 2 1 distribution.
We recognize this method does not include the variance component and is thus, not the same model that we are proposing.The likelihood ratio test for GLM described above provides an idea about the relative important of various fixed effects covariates that may be influential for light, moderate or heavy precipitation.The above procedure may be supplemented by a multiple testing correction procedure, if needed.The details of this analysis is available from the authors.Also note that in this part of the analysis we did not include random effects, owing to a lack of viable and theoretically justifiable testing procedure when a random effect is present.Inclusion of random effects are likely to reduce variance attributed to noise, thus typically increasing significance levels.
We chose the GLM with all fixed covariates to provide a test of goodness-of-fit based on residual deviance being asymptotically χ 2 .For details, refer to Faraway (2006).This compares the fitted model to the saturated model which contains one parameter for each observation.Failure to reject the null hypothesis of this test indicates a lack-of-fit.
Finally, a test of the variance component in the GLMM fit by glmer is done using a LRT with a nonstandard asymptotic distribution.Because our models have a single variance component, the asymptotic distribution for the LRT corresponds half of the p-value obtained from the χ 2 1 distribution as noted by Zhang and Lin (2008).We only do this test for glmer since the other two methods do not use maximum likelihood, although it should be noted the likelihood produced by glmer is only approximate.Overall, we take a cautious view on the interpretation of these tests.

Discussion of Rainfall Models
Results of significance testing for each of the models can be found in Table 6.All covariates were significant in the light rainfall model in the majority of the years.However, 38% of the years showed lack-of-fit based on the deviance test.The moderate and extreme rainfall models showed no lack-of-fit, but had far fewer significant covariates over the years.
Clearly, maximum daily temperature was important in all three levels of rainfall aligning with the Clausius-Clapeyron equation regarding water vapor capacity of the atmosphere.Minimum temperature was significant in most years for light and moderate rainfall, but was only significant in a minority 595 of the extreme rainfall models.
Elevation was also significant in many years for all rainfall levels.This aligns with the physical explanations of warm moist air cooling at higher altitudes to produce precipitation.
Latitude and longitude were both significant in most light 600 rainfall years.Moderate and extreme rainfall did not indicate latitude as significant in most years.Longitude was significant in just over half of the extreme models.Coefficient estimates for latitude indicated the probabilities of rain increasing going south to north.Longitude estimates were mostly 605 negative indicating a decreasing probability of rainfall going west to east.DMI was significant for the majority of light rainfall models; however, it was significant in very few of the moderate and extreme models.This corresponds to the DMI influence 610 in non-ENSO drought years as hypothesized in previous literature.
∆T T was significant in most years for all three rainfall levels corroborating the hypothesis that it is instrumental in monsoon circulation.

615
The Nino 3.4 anomaly index was significant in the majority of light rainfall models, but in less than 20% of both moderate and extreme models.This could be related to the possible weakening of the relationship between ENSO and the Indian monsoon as noted in Chang et al. (2001) but may 620 also be a function of the other covariates included in the modeling.
The station variance component was significant in nearly all of the light rainfall models.One can note from Fig. 2a, there is less variability in general in light rainfall even though 625 most years are significant.The variance was significant in about half of the moderate and a quarter of the extreme rainfall models.As seen in Figs.2b-2c, these models tended to have higher variability than light rainfall even though fewer years were significant.The variance component does provide 630 additional explanation for the rainfall variability and thus, vets the methodology use in this application.This verifies the thesis of this paper: that a significant portion of the variability in any precipitation category is a random component that is distinguishable from random noise variability.

Estimation Method Performance
The coefficient estimates over the time period for all fixed effects at each level of rainfall are depicted in Figs.3-5.The three estimation methods tended to produce different answers on at least a few of the coefficients in each of the 640 models.The best agreement amongst all methods occurred in the estimates for maximum temperature (moderate, extreme), longitude (light, extreme), and the variance components (all).
Estimates for glmer and MSIM tended to agree more often than either agreed with dclone.However, they were fairly different in magnitude for several covariates and did not always trend with each other.Light rainfall models displayed slightly more disagreement more than moderate and extreme rainfall.
Standard deviations from dclone estimates indicated that the algorithm had likely not converged for all parameters in the 10000 samples taken from the posterior.As mentioned in §3.3, one of the issues with using this method is difficulty in assessing convergence.It would likely require a much larger sample to provide suitable answers in this application.Based on this, we would say the dclone results were mostly inconclusive in this application.
The outcome of this application indicated glmer and MSIM provided more reasonable estimates, however, a longer run of dclone may also be useful.The three methods are representative of distinct statistical paradigms of estimation including approximate likelihood, Bayesian, and method of moments.Each of the methods uses a different algorithm and different assumptions, thus, we recommend use of multiple methods when applying GLMM in an application.

Conclusions
Outcomes for Indian monsoon rainfall coincide with results in the Indian monsoon literature providing evidence of the usefulness of GLMM.Physical constructs are preserved by the models demonstrated by the importance of elevation, maximum temperature, and ∆T T in all levels of rainfall.Random effects are significant in several of the models indicating promise of modeling some of the unobservable and complicated interactions that underly climate patterns.
The GLMM methods explored in this article purposely included several styles of estimation including approximate likelihood, Bayesian, and method of moments type estimators.Each exhibited some drawbacks, thus, use of at least two out of three of the best methods, glmer, MSIM, and dclone, in the context of any application will help verify consistency of estimates.Use of multiple methods will provide researchers with higher confidence in results and will be more robust to limitations of any of the individual methods.
Since the relevance of GLMM in this context has been established, climate model output, such as that of CMIP5, will be explored to gain deeper intuition of the nature of this random effect.Further work on GLMMs may include studying other proposed drivers of Indian monsoons in their contributions to fixed or random effects.Additional random effects that include spatial correlation will be tested in future work.
We also note that this model could be pursued in the future as a multinomial logit model.
It was suggested that Normalized Difference Vegetation Index (NDVI) may be a useful covariate for understanding feedback of vegetation on precipitation.However, data cov-erage was limited, thus, it was not included in our models.It will be examined more closely in future modeling efforts.
Providing improvements to the GLMM estimation methods is another open research area.One limitation of GLMM, as presented here, is the reliance of modeling random ef-700 fects as normal.Expanding the possible distributions of random effects to include extreme value distributions would be a breakthrough in mixed modeling.Providing more conclusive tests of significance for fixed and random effects is also important.

Table 5 .
Total System Time (in seconds) Results for nokomis

Table 6 .
This table indicates the percentage of significant p-values at α = 0.05 level for each of the 1973-2013 models.p-values for fixed coefficients and goodness-of-fit test are from LRTs on GLM fits.p-values for the variance components are from LRTs that compare GLM and glmer fits.