Source code for naginterfaces.library.examples.correg.lmm_init_combine_ex

#!/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 )