MixedModels.GHnorm — FunctionGHnorm(k::Int)Return the (unique) GaussHermiteNormalized{k} object.
The function values are stored (memoized) when first evaluated. Subsequent evaluations for the same k have very low overhead.
gausshermitenorm (generic function with 1 method)providing
gausshermitenorm(3)([-1.7320508075688739, 1.1102230246251565e-15, 1.7320508075688774], [0.16666666666666743, 0.6666666666666657, 0.16666666666666677])The weights and positions are often shown as a lollipop plot. For the 9th order rule these are
gh9=gausshermitenorm(9)
-plot(x=gh9[1], y=gh9[2], Geom.hair, Geom.point, Guide.ylabel("Weight"), Guide.xlabel(""))
Notice that the magnitudes of the weights drop quite dramatically away from zero, even on a logarithmic scale
plot(
+plot(x=gh9[1], y=gh9[2], Geom.hair, Geom.point, Guide.ylabel("Weight"), Guide.xlabel(""))
Notice that the magnitudes of the weights drop quite dramatically away from zero, even on a logarithmic scale
plot(
x=gh9[1], y=gh9[2], Geom.hair, Geom.point,
Scale.y_log2, Guide.ylabel("Weight (log scale)"),
Guide.xlabel(""),
-)
The definition of MixedModels.GHnorm is similar to the gausshermitenorm function with some extra provisions for ensuring symmetry of the abscissae and the weights and for caching values once they have been calculated.
using MixedModels
+)
The definition of MixedModels.GHnorm is similar to the gausshermitenorm function with some extra provisions for ensuring symmetry of the abscissae and the weights and for caching values once they have been calculated.
MixedModels.GHnorm — FunctionGHnorm(k::Int)Return the (unique) GaussHermiteNormalized{k} object.
The function values are stored (memoized) when first evaluated. Subsequent evaluations for the same k have very low overhead.
using MixedModels
GHnorm(3)MixedModels.GaussHermiteNormalized{3}([-1.7320508075688772, 0.0, 1.7320508075688772], [0.16666666666666666, 0.6666666666666666, 0.16666666666666666])By the properties of the normal distribution, when $\mathcal{X}\sim\mathscr{N}(\mu, \sigma^2)$
\[\mathbb{E}[g(x)] \approx \sum_{i=1}^k g(\mu + \sigma z_i)\,w_i\]
For example, $\mathbb{E}[\mathcal{X}^2]$ where $\mathcal{X}\sim\mathcal{N}(2, 3^2)$ is
μ = 2; σ = 3; ghn3 = GHnorm(3);
sum(@. ghn3.w * abs2(μ + σ * ghn3.z)) # should be μ² + σ² = 1313.0(In general a dot, '.', after the function name in a function call, as in abs2.(...), or before an operator creates a fused vectorized evaluation in Julia. The macro @. has the effect of vectorizing all operations in the subsequent expression.)
Application to a model for contraception use
A binary response is a "Yes"/"No" type of answer. For example, in a 1989 fertility survey of women in Bangladesh (reported in Huq, N. M. and Cleland, J., 1990) one response of interest was whether the woman used artificial contraception. Several covariates were recorded including the woman's age (centered at the mean), the number of live children the woman has had (in 4 categories: 0, 1, 2, and 3 or more), whether she lived in an urban setting, and the district in which she lived. The version of the data used here is that used in review of multilevel modeling software conducted by the Center for Multilevel Modelling, currently at University of Bristol (http://www.bristol.ac.uk/cmm/learning/mmsoftware/data-rev.html). These data are available as the :contra dataset.
contra = DataFrame(MixedModels.dataset(:contra))
describe(contra)5×7 DataFrame
| Row | variable | mean | min | median | max | nmissing | eltype |
|---|---|---|---|---|---|---|---|
| Symbol | Union… | Any | Union… | Any | Int64 | DataType | |
| 1 | dist | D01 | D61 | 0 | String | ||
| 2 | urban | N | Y | 0 | String | ||
| 3 | livch | 0 | 3+ | 0 | String | ||
| 4 | age | 0.00204757 | -13.56 | -1.56 | 19.44 | 0 | Float64 |
| 5 | use | N | Y | 0 | String |
A smoothed scatterplot of contraception use versus age
plot(contra, x=:age, y=:use, Geom.smooth, Guide.xlabel("Centered age (yr)"),
- Guide.ylabel("Contraception use"))
shows that the proportion of women using artificial contraception is approximately quadratic in age.
A model with fixed-effects for age, age squared, number of live children and urban location and with random effects for district, is fit as
const form1 = @formula use ~ 1 + age + abs2(age) + livch + urban + (1|dist);
+ Guide.ylabel("Contraception use"))
shows that the proportion of women using artificial contraception is approximately quadratic in age.
A model with fixed-effects for age, age squared, number of live children and urban location and with random effects for district, is fit as
const form1 = @formula use ~ 1 + age + abs2(age) + livch + urban + (1|dist);
m1 = fit(MixedModel, form1, contra, Bernoulli(), fast=true)Generalized Linear Mixed Model fit by maximum likelihood (nAGQ = 1)
use ~ 1 + age + :(abs2(age)) + livch + urban + (1 | dist)
Distribution: Bernoulli{Float64}
@@ -86,15 +86,15 @@
-
+
+
u₁
-
+
-5
@@ -111,52 +111,52 @@
-
-
-
-
+
+
+
+
-
-
+
+
-
+
-
+
-
+
-
+
-
+
-
-
+
+
-
+
-
+
-
+
-
-
-
+
+
+
@@ -164,40 +164,40 @@
-
-
+
+
0
-
+
100
-
+
200
-
+
300
-
+
400
-
+
500
-
-
+
+
Deviance contribution
@@ -205,7 +205,7 @@
-
+
@@ -227,15 +227,15 @@
-
+
+
u₃
-
+
-5
@@ -252,49 +252,49 @@
-
-
-
-
+
+
+
+
-
-
+
+
-
+
-
+
-
+
-
+
-
-
+
+
-
+