#!/usr/bin/env python3
"``naginterfaces.library.correg.lmm_init_combine`` Python Example."
# NAG Copyright 2017-2019.
# pylint: disable=invalid-name,too-many-locals,too-many-statements
from collections import namedtuple
import warnings
import numpy as np
from naginterfaces.base import utils
from naginterfaces.library import blgm, correg
[docs]def main():
"""
Example for :func:`naginterfaces.library.correg.lmm_init_combine`
including calls to
:func:`naginterfaces.library.blgm.lm_submodel`,
:func:`naginterfaces.library.blgm.lm_describe_data`,
:func:`naginterfaces.library.correg.lmm_init`,
:func:`naginterfaces.library.correg.lmm_fit` and
:func:`naginterfaces.library.blgm.handle_free`
Multi-level linear mixed effects regression model using
restricted maximum likelihood.
The data used in this example was simulated.
>>> main()
naginterfaces.library.correg.lmm_init_combine Python Example Results.
Linear mixed effects regression model using REML.
Random Parameter Estimates
==========================
Estimate Standard Label
Error
0.683 0.506 V7 . V12 (Lvl 1)
...
0.504 2.693 V4 (Lvl 3) . V11 (Lvl 3) . V10 (Lvl 2) . V12 (Lvl 3)
Fixed Parameter Estimates
=========================
Estimate Standard Label
Error
1.643 2.460 Intercept
-1.622 0.855 V1 (Lvl 2)
-2.482 1.142 V2 (Lvl 2)
0.462 1.214 V2 (Lvl 3)
Variance Components
===================
Estimate Label
0.563 V7 . V12
5.820 V8 . V12
10.860 V9 . V12
19.628 V5 . V11 . V12
40.534 V6 . V11 . V12
36.323 V3 . V11 . V10 . V12
12.451 V4 . V11 . V10 . V12
Sigma^2 = 0.003
-2 Log Likelihood = 608.195
"""
def dataset_meta_data():
"""
Return some meta data for the dataset
"""
# The levels and variable names associated with the independent
# variables:
levels = [2, 3, 2, 3, 2, 3, 1, 4, 5, 2, 3, 3]
vnames = ['V' + str(x + 1) for x in range(len(levels))]
# Number of observations (in this context, this is not really used, so
# can take any value) we have set it here to the number of observations
# in the full dataset because we readily know it
nobs = 90
return namedtuple(
'meta_data',
['nobs', 'levels', 'vnames'])(
nobs, levels, vnames
)
def return_dataset(fobs=None, lobs=None):
"""
Return data points for observations fobs:lobs
"""
# The dependent data:
y = [3.1100, 2.8226, 7.4543, 4.4313, 6.1543, -0.1783, 4.6748,
7.0667, 1.4262, 7.7290, -2.1806, 6.8419, 1.2590, 8.8405,
6.1657, -4.5605, -1.2367, -12.2932, -2.3374, 0.0716,
0.1895, 1.5608, -0.8529, -4.1169, 3.9977, -8.1277,
-4.9656, -0.6428, -5.5152, -5.5657, 14.8177, 16.9783,
13.8966, 14.8166, 19.3640, 9.5299, 12.0102, 6.1551,
-1.7048, 2.7640, 2.8065, 0.0974, -7.8080, -18.0450,
-2.8199, 8.9893, 3.7978, -6.3493, 8.1411, -7.5483,
-0.4600, -3.2135, -6.6562, 5.1267, 3.5592, -4.4420,
-8.5965, -6.3187, -7.8953, -10.1383, -7.8850, 23.2001,
5.5829, -4.3698, 2.1274, -2.7184, -17.9128, -1.2708,
-24.2735, -14.7374, 0.1713, 8.0006, 1.2100, 3.3307,
-22.6713, 7.5562, -7.0694, 3.7159, -4.3135, -14.5577,
-12.5107, 4.7708, 13.2797, -6.3243, -7.0549, -9.2713,
-18.7788, -7.7230, -22.7230, -11.6609]
# The independent data:
dat = [
[1.0, 3.0, 2.0, 1.0, 2.0, 2.0, -0.3160, 4.0, 2.0, 1.0, 1.0, 1.0],
[1.0, 1.0, 1.0, 3.0, 1.0, 2.0, -1.3377, 1.0, 4.0, 1.0, 1.0, 1.0],
[1.0, 3.0, 1.0, 3.0, 1.0, 3.0, -0.7610, 4.0, 2.0, 1.0, 1.0, 1.0],
[2.0, 3.0, 2.0, 1.0, 1.0, 3.0, -2.2976, 4.0, 2.0, 1.0, 1.0, 1.0],
[2.0, 2.0, 1.0, 3.0, 2.0, 3.0, -0.4263, 2.0, 1.0, 1.0, 1.0, 1.0],
[2.0, 1.0, 2.0, 3.0, 1.0, 3.0, 1.4067, 3.0, 3.0, 2.0, 1.0, 1.0],
[2.0, 3.0, 2.0, 1.0, 2.0, 1.0, -1.4669, 1.0, 2.0, 2.0, 1.0, 1.0],
[1.0, 1.0, 1.0, 3.0, 2.0, 3.0, 0.4717, 2.0, 4.0, 2.0, 1.0, 1.0],
[1.0, 3.0, 2.0, 3.0, 2.0, 1.0, 0.4436, 1.0, 3.0, 2.0, 1.0, 1.0],
[1.0, 1.0, 1.0, 2.0, 2.0, 3.0, -0.5950, 3.0, 4.0, 2.0, 1.0, 1.0],
[1.0, 3.0, 1.0, 3.0, 1.0, 1.0, -1.7981, 4.0, 2.0, 1.0, 2.0, 1.0],
[2.0, 3.0, 1.0, 2.0, 1.0, 1.0, 0.2397, 1.0, 4.0, 1.0, 2.0, 1.0],
[1.0, 2.0, 2.0, 1.0, 2.0, 3.0, 0.4742, 1.0, 1.0, 1.0, 2.0, 1.0],
[2.0, 2.0, 2.0, 2.0, 2.0, 3.0, 0.6888, 3.0, 1.0, 1.0, 2.0, 1.0],
[2.0, 1.0, 2.0, 3.0, 1.0, 3.0, -1.0616, 3.0, 5.0, 1.0, 2.0, 1.0],
[1.0, 2.0, 2.0, 2.0, 2.0, 1.0, -0.5356, 1.0, 3.0, 2.0, 2.0, 1.0],
[1.0, 3.0, 2.0, 2.0, 1.0, 1.0, -1.2963, 2.0, 5.0, 2.0, 2.0, 1.0],
[1.0, 2.0, 2.0, 1.0, 2.0, 2.0, -1.5389, 3.0, 2.0, 2.0, 2.0, 1.0],
[2.0, 3.0, 1.0, 1.0, 2.0, 2.0, -0.6408, 2.0, 1.0, 2.0, 2.0, 1.0],
[1.0, 2.0, 2.0, 2.0, 1.0, 1.0, 0.6574, 1.0, 1.0, 2.0, 2.0, 1.0],
[2.0, 1.0, 1.0, 1.0, 1.0, 3.0, 0.9259, 1.0, 2.0, 1.0, 3.0, 1.0],
[2.0, 2.0, 2.0, 1.0, 2.0, 2.0, 1.5080, 3.0, 1.0, 1.0, 3.0, 1.0],
[2.0, 3.0, 1.0, 1.0, 1.0, 3.0, 2.5821, 2.0, 3.0, 1.0, 3.0, 1.0],
[1.0, 2.0, 2.0, 1.0, 2.0, 3.0, 0.4102, 1.0, 4.0, 1.0, 3.0, 1.0],
[2.0, 1.0, 2.0, 3.0, 2.0, 2.0, 0.7839, 2.0, 5.0, 1.0, 3.0, 1.0],
[1.0, 2.0, 2.0, 3.0, 2.0, 1.0, -1.8812, 4.0, 2.0, 2.0, 3.0, 1.0],
[1.0, 2.0, 1.0, 3.0, 2.0, 3.0, 0.7770, 4.0, 1.0, 2.0, 3.0, 1.0],
[2.0, 2.0, 1.0, 2.0, 1.0, 3.0, 0.2590, 3.0, 1.0, 2.0, 3.0, 1.0],
[2.0, 3.0, 2.0, 2.0, 2.0, 3.0, -0.9250, 3.0, 3.0, 2.0, 3.0, 1.0],
[2.0, 2.0, 1.0, 3.0, 2.0, 3.0, -0.4831, 1.0, 5.0, 2.0, 3.0, 1.0],
[2.0, 2.0, 1.0, 3.0, 1.0, 3.0, 0.5046, 3.0, 3.0, 1.0, 1.0, 2.0],
[2.0, 1.0, 1.0, 2.0, 2.0, 1.0, -0.6903, 2.0, 1.0, 1.0, 1.0, 2.0],
[1.0, 3.0, 2.0, 2.0, 2.0, 1.0, 1.6166, 2.0, 5.0, 1.0, 1.0, 2.0],
[2.0, 2.0, 2.0, 2.0, 1.0, 3.0, 0.2778, 2.0, 3.0, 1.0, 1.0, 2.0],
[2.0, 3.0, 2.0, 2.0, 1.0, 2.0, 1.9586, 4.0, 2.0, 1.0, 1.0, 2.0],
[1.0, 3.0, 1.0, 1.0, 1.0, 3.0, 1.0506, 2.0, 5.0, 2.0, 1.0, 2.0],
[2.0, 1.0, 1.0, 3.0, 2.0, 3.0, 0.4871, 1.0, 1.0, 2.0, 1.0, 2.0],
[2.0, 1.0, 2.0, 3.0, 2.0, 1.0, 2.0891, 4.0, 4.0, 2.0, 1.0, 2.0],
[1.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.4338, 4.0, 3.0, 2.0, 1.0, 2.0],
[1.0, 1.0, 2.0, 3.0, 1.0, 2.0, -1.1196, 3.0, 4.0, 2.0, 1.0, 2.0],
[1.0, 3.0, 1.0, 1.0, 2.0, 1.0, 0.3367, 3.0, 2.0, 1.0, 2.0, 2.0],
[2.0, 2.0, 1.0, 3.0, 1.0, 1.0, 0.1092, 2.0, 2.0, 1.0, 2.0, 2.0],
[1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 0.4007, 4.0, 1.0, 1.0, 2.0, 2.0],
[2.0, 3.0, 1.0, 1.0, 1.0, 2.0, 0.1460, 3.0, 5.0, 1.0, 2.0, 2.0],
[2.0, 1.0, 2.0, 3.0, 1.0, 3.0, -0.3877, 3.0, 4.0, 1.0, 2.0, 2.0],
[1.0, 1.0, 1.0, 2.0, 2.0, 1.0, 0.6957, 4.0, 3.0, 2.0, 2.0, 2.0],
[2.0, 1.0, 1.0, 1.0, 2.0, 1.0, -0.4664, 3.0, 3.0, 2.0, 2.0, 2.0],
[1.0, 1.0, 1.0, 1.0, 2.0, 3.0, 0.2067, 2.0, 4.0, 2.0, 2.0, 2.0],
[2.0, 1.0, 2.0, 1.0, 1.0, 2.0, 0.4112, 1.0, 4.0, 2.0, 2.0, 2.0],
[2.0, 2.0, 1.0, 1.0, 1.0, 2.0, -1.3734, 3.0, 3.0, 2.0, 2.0, 2.0],
[2.0, 1.0, 2.0, 3.0, 1.0, 3.0, 0.7065, 1.0, 3.0, 1.0, 3.0, 2.0],
[1.0, 2.0, 2.0, 2.0, 1.0, 2.0, 1.3628, 4.0, 2.0, 1.0, 3.0, 2.0],
[2.0, 1.0, 2.0, 2.0, 2.0, 3.0, -0.5052, 4.0, 5.0, 1.0, 3.0, 2.0],
[2.0, 1.0, 1.0, 1.0, 2.0, 1.0, -1.3457, 2.0, 5.0, 1.0, 3.0, 2.0],
[1.0, 1.0, 2.0, 1.0, 2.0, 3.0, -1.8022, 3.0, 4.0, 1.0, 3.0, 2.0],
[2.0, 3.0, 1.0, 2.0, 1.0, 1.0, 0.0116, 2.0, 4.0, 2.0, 3.0, 2.0],
[2.0, 2.0, 1.0, 3.0, 2.0, 3.0, -0.9075, 1.0, 3.0, 2.0, 3.0, 2.0],
[2.0, 2.0, 2.0, 2.0, 2.0, 3.0, -1.4707, 1.0, 1.0, 2.0, 3.0, 2.0],
[2.0, 2.0, 1.0, 1.0, 2.0, 1.0, -1.2938, 2.0, 3.0, 2.0, 3.0, 2.0],
[1.0, 3.0, 1.0, 3.0, 2.0, 2.0, -1.1660, 4.0, 4.0, 2.0, 3.0, 2.0],
[1.0, 2.0, 1.0, 1.0, 2.0, 3.0, 0.0397, 4.0, 4.0, 1.0, 1.0, 3.0],
[1.0, 3.0, 1.0, 2.0, 1.0, 3.0, -0.5987, 3.0, 2.0, 1.0, 1.0, 3.0],
[2.0, 3.0, 2.0, 2.0, 1.0, 1.0, 0.6683, 3.0, 3.0, 1.0, 1.0, 3.0],
[2.0, 2.0, 1.0, 1.0, 2.0, 2.0, -0.0106, 1.0, 3.0, 1.0, 1.0, 3.0],
[1.0, 2.0, 1.0, 3.0, 2.0, 2.0, 0.5885, 1.0, 3.0, 1.0, 1.0, 3.0],
[1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 0.4555, 1.0, 5.0, 2.0, 1.0, 3.0],
[2.0, 2.0, 2.0, 1.0, 1.0, 2.0, 0.6502, 4.0, 3.0, 2.0, 1.0, 3.0],
[1.0, 1.0, 1.0, 3.0, 1.0, 1.0, -0.1601, 1.0, 3.0, 2.0, 1.0, 3.0],
[2.0, 2.0, 1.0, 3.0, 2.0, 3.0, 1.6910, 1.0, 1.0, 2.0, 1.0, 3.0],
[2.0, 2.0, 2.0, 3.0, 1.0, 2.0, 0.1053, 4.0, 4.0, 2.0, 1.0, 3.0],
[2.0, 1.0, 2.0, 3.0, 2.0, 2.0, -0.4037, 3.0, 4.0, 1.0, 2.0, 3.0],
[1.0, 3.0, 2.0, 3.0, 1.0, 3.0, -0.5853, 3.0, 2.0, 1.0, 2.0, 3.0],
[2.0, 3.0, 2.0, 1.0, 1.0, 1.0, -0.3037, 1.0, 3.0, 1.0, 2.0, 3.0],
[1.0, 3.0, 1.0, 1.0, 2.0, 2.0, -0.0774, 1.0, 4.0, 1.0, 2.0, 3.0],
[2.0, 3.0, 1.0, 2.0, 2.0, 1.0, 0.4733, 4.0, 5.0, 1.0, 2.0, 3.0],
[1.0, 3.0, 2.0, 2.0, 1.0, 2.0, -0.0354, 4.0, 2.0, 2.0, 2.0, 3.0],
[1.0, 3.0, 2.0, 2.0, 1.0, 1.0, -0.6640, 2.0, 1.0, 2.0, 2.0, 3.0],
[2.0, 3.0, 1.0, 3.0, 1.0, 1.0, 0.0335, 4.0, 4.0, 2.0, 2.0, 3.0],
[1.0, 2.0, 2.0, 2.0, 1.0, 3.0, 0.1351, 1.0, 1.0, 2.0, 2.0, 3.0],
[1.0, 1.0, 2.0, 1.0, 2.0, 3.0, -0.5951, 3.0, 4.0, 2.0, 2.0, 3.0],
[2.0, 2.0, 2.0, 3.0, 1.0, 3.0, 0.2735, 3.0, 2.0, 1.0, 3.0, 3.0],
[2.0, 2.0, 1.0, 1.0, 1.0, 3.0, 0.3157, 1.0, 2.0, 1.0, 3.0, 3.0],
[2.0, 2.0, 2.0, 1.0, 1.0, 1.0, -1.0843, 2.0, 3.0, 1.0, 3.0, 3.0],
[1.0, 2.0, 2.0, 1.0, 2.0, 2.0, -0.0836, 4.0, 2.0, 1.0, 3.0, 3.0],
[2.0, 1.0, 2.0, 1.0, 1.0, 2.0, -0.2884, 2.0, 1.0, 1.0, 3.0, 3.0],
[2.0, 3.0, 2.0, 3.0, 2.0, 3.0, -0.1006, 1.0, 2.0, 2.0, 3.0, 3.0],
[1.0, 3.0, 1.0, 2.0, 2.0, 3.0, 0.5710, 1.0, 3.0, 2.0, 3.0, 3.0],
[1.0, 1.0, 2.0, 1.0, 1.0, 2.0, 0.2776, 2.0, 3.0, 2.0, 3.0, 3.0],
[2.0, 3.0, 2.0, 2.0, 1.0, 3.0, -0.7561, 4.0, 4.0, 2.0, 3.0, 3.0],
[1.0, 2.0, 2.0, 2.0, 1.0, 2.0, 1.5549, 1.0, 4.0, 2.0, 3.0, 3.0],
]
tfobs = 0 if fobs is None else fobs
tlobs = len(y) if lobs is None else lobs
return namedtuple(
'data', ['dat', 'y'],)(
dat[tfobs:tlobs], y[tfobs:tlobs]
)
def get_vinfo(what, hlmm):
"""
Call blgm.lm_submodel to get descriptive information on
estimates returned from fitting the linear mixed effects
regression model defined in hlmm
"""
# supply a dummy submodel handle
hform = utils.Handle()
# Find how long vinfo needs to be
lvinfo = 3
with warnings.catch_warnings():
warnings.simplefilter('ignore', utils.NagAlgorithmicWarning)
subres = blgm.lm_submodel(what, hform, hlmm, 0, 0, lvinfo)
lvinfo = subres.vinfo[2]
# Get and return vinfo
subres = blgm.lm_submodel(what, hform, hlmm, 0, 0, lvinfo)
return subres.vinfo
def get_next_label(vinfo, vnames=None, fixed_len=None):
"""
Return a label based on the information supplied in vinfo
Routine should be called for all estimates described by
vinfo, in order. vinfo is updated in place in preparation
for the next call
"""
k = 0
b = vinfo[k]
k += 1
if b < 0:
# the corresponding column of the design matrix was padding
label = ''
elif b == 0:
label = 'Intercept'
else:
label = ''
for _i in range(b):
# Extract the information from VINFO
# (we ignore the contrast information)
vi, li, _ = vinfo[k:k+3]
k += 3
# Variable name
if vnames is None:
this_vname = 'V' + str(vi)
else:
# vi'th variable supplied when describing the dataset
this_vname = vnames[vi - 1]
# Construct the label
if label != '':
label += ' . '
label += this_vname
if li > 0:
label += ' (Lvl ' + str(li) + ')'
if fixed_len is not None:
label = label.ljust(fixed_len)
# shift vinfo in place, so it is ready to produce the
# next label
np.copyto(vinfo, np.roll(vinfo, -k))
return label
print(
'naginterfaces.library.correg.lmm_init_combine Python Example Results.'
)
print('Linear mixed effects regression model using REML.')
# Get the meta data for the dataset
meta = dataset_meta_data()
hddesc = utils.Handle()
blgm.lm_describe_data(hddesc, meta.nobs, meta.levels, vnames=meta.vnames)
# Specify the fixed part of the model:
hfixed = utils.Handle()
blgm.lm_formula(hfixed, 'V1 + V2')
# Specify the random part of the model:
hrndm = [utils.Handle() for i in range(3)]
blgm.lm_formula(hrndm[0], '-1 + V7 + V8 + V9 / SUBJECT = V12')
blgm.lm_formula(hrndm[1], '-1 + V5 + V6 / SUBJECT = V12.V11')
blgm.lm_formula(hrndm[2], '-1 + V3 + V4 / SUBJECT = V12.V11.V10')
# Handle to hold a description of the mixed effects model
hlmm = utils.Handle()
# In this contrived example we will be processing the data
# in various sized blocks
blk_size = [10, 31, 40, 9]
# Process all the blocks of data
fobs = 0
for nobs in blk_size:
lobs = fobs + nobs
# Get the data for observations fobs:lobs
this_blk = return_dataset(fobs, lobs)
# Process the current block of data
pdat = correg.lmm_init(hlmm, hddesc, hfixed, this_blk.y,
this_blk.dat, hrndm=hrndm)
if fobs == 0:
xcomm = pdat.comm.copy()
else:
correg.lmm_init_combine(hlmm, xcomm, pdat.comm)
fobs = lobs
# Initial values for the variance components
# the number of variance components (nvpr) is not data
# dependent, so we can get it from pdat for any of the data blocks
gamma = [-1.0] * (pdat.nvpr + 1)
# Fit the model
fit = correg.lmm_fit(hlmm, gamma, xcomm)
# Get descriptive information on the estimates for the
# fixed parameters, random parameters and variance components
fixed_vinfo = get_vinfo('X', hlmm)
random_vinfo = get_vinfo('Z', hlmm)
vc_vinfo = get_vinfo('V', hlmm)
# The number of random and fixed factors for the full dataset
# (the ones returned by lmm_init are for a particular subset of data):
nff = correg.optget('NFF', xcomm)['value']
nrf = correg.optget('NRF', xcomm)['value']
rnlsv = correg.optget('RNLSV', xcomm)['value']
fnlsv = correg.optget('FNLSV', xcomm)['value']
nxx = nff*fnlsv
nzz = nrf*rnlsv
# Display the results
pad = ' '
print('Random Parameter Estimates')
print('==========================')
print(' Estimate Standard Label')
print(' Error')
for p in range(nzz):
label = get_next_label(random_vinfo, meta.vnames)
if label != '':
print('{:8.3f}'.format(fit.b[p]) + '{:8.3f}'.format(fit.se[p]) +
pad + label)
print('Fixed Parameter Estimates')
print('=========================')
print(' Estimate Standard Label')
print(' Error')
for p in range(nxx):
label = get_next_label(fixed_vinfo, meta.vnames)
if label != '':
print('{:8.3f}'.format(fit.b[p+nzz]) +
'{:8.3f}'.format(fit.se[p+nzz]) +
pad + label)
print('Variance Components')
print('===================')
print(' Estimate Label')
for p in range(pdat.nvpr):
print('{:9.3f}'.format(fit.gamma[p]) + pad +
get_next_label(vc_vinfo, meta.vnames))
print('Sigma^2 = {:15.3f}'.format(fit.gamma[pdat.nvpr]))
print('-2 Log Likelihood = {:15.3f}'.format(fit.lnlike))
# Free the rest of the blgm handles:
map(blgm.handle_free, [hfixed, hlmm] + hrndm)
if __name__ == '__main__':
import doctest
import sys
sys.exit(
doctest.testmod(
None, verbose=True, report=False,
optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
).failed
)
