From 644a66d63e2dd00cffc0be31953a03530136020d Mon Sep 17 00:00:00 2001 From: Charles Driver Date: Tue, 28 May 2024 10:08:48 +0200 Subject: [PATCH] .. --- DESCRIPTION | 2 +- R/fitIRT.R | 107 ++++++++++------ bigIRT.Rproj | 1 - inst/stan/{irtmml.stan => irtmml.bak} | 0 readme.Rmd | 84 ++++++++++--- readme.md | 119 +++++++++++++----- readme_files/figure-gfm/unnamed-chunk-1-1.png | Bin 7934 -> 11368 bytes readme_files/figure-gfm/unnamed-chunk-1-2.png | Bin 0 -> 11696 bytes readme_files/figure-gfm/unnamed-chunk-1-3.png | Bin 0 -> 11838 bytes 9 files changed, 223 insertions(+), 90 deletions(-) rename inst/stan/{irtmml.stan => irtmml.bak} (100%) create mode 100644 readme_files/figure-gfm/unnamed-chunk-1-2.png create mode 100644 readme_files/figure-gfm/unnamed-chunk-1-3.png diff --git a/DESCRIPTION b/DESCRIPTION index b9499b6..f4b7508 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: bigIRT Title: Fits item response theory models to big data -Version: 0.1.6 +Version: 0.1.7 Authors@R: person(given = "Charles", family = "Driver", diff --git a/R/fitIRT.R b/R/fitIRT.R index 12f7d77..ecaebd5 100644 --- a/R/fitIRT.R +++ b/R/fitIRT.R @@ -300,44 +300,67 @@ dropPerfectScores <- function(dat,scoreref.='score',itemref.='Item',idref.='id', # } -#' Title +#' Fit a binary Item Response Theory (IRT) model +#' +#' This function fits a binary Item Response Theory (IRT) model using various parameters and options. +#' +#' @param dat A data frame containing the data to be analyzed. +#' @param score Character. The name of the column in \code{dat} representing the response. Default is 'score'. +#' @param id Character. The name of the column in \code{dat} representing the individual. Default is 'id'. +#' @param item Character. The name of the column in \code{dat} representing the item. Default is 'Item'. +#' @param scale Character. The name of the column in \code{dat} representing the scale of each item. Default is 'Scale'. +#' @param pl Integer. The number of parameters for the logistic model (1PL, 2PL, 3PL, or 4PL). Default is 1. +#' @param personDat Data frame. Any fixed ability data for persons. Default is NA. +#' @param personPreds Character vector. Names of predictors for person parameters found in data. Default is an empty character vector. +#' @param itemDat Data frame. Any fixed item data. Default is NA. +#' @param AitemPreds Character vector. Names of predictors for item discrimination parameters. Default is an empty character vector. +#' @param BitemPreds Character vector. Names of predictors for item difficulty parameters. Default is an empty character vector. +#' @param CitemPreds Character vector. Names of predictors for item guessing parameters. Default is an empty character vector. +#' @param DitemPreds Character vector. Names of predictors for item upper asymptote (slipping) parameters. Default is an empty character vector. +#' @param itemSpecificBetas Logical. Whether to allow item-specific betas for covariate effects, or simply estimate one effect per covariate. Default is FALSE. +#' @param betaScale Numeric. Scale of the prior for beta parameters. Default is 10. +#' @param invspAMeandat Numeric. Mean for the prior distribution of the raw discrimination parameters, +#' which subsequently have a 'softplus' log(1+exp(x)) applied. Default is 0.542, giving a mean for A pars of ~ 1. +#' @param invspASD Numeric. Standard deviation for the prior distribution of the raw discrimination parameters. Default is 2. +#' @param BMeandat Numeric. Mean for the prior distribution of the item difficulty parameters. Default is 0. +#' @param BSD Numeric. Standard deviation for the prior distribution of the item difficulty parameters. Default is 10. +#' @param logitCMeandat Numeric. Mean for the prior distribution of the item guessing parameters (on logit scale). Default is -4. +#' @param logitCSD Numeric. Standard deviation for the prior distribution of the item guessing parameters (on logit scale). Default is 2. +#' @param logitDMeandat Numeric. Mean for the prior distribution of the item upper asymptote parameters (on logit scale). Default is 4. +#' @param logitDSD Numeric. Standard deviation for the prior distribution of the item upper asymptote parameters (on logit scale). Default is 2. +#' @param AbilityMeandat Numeric array. Mean for the prior distribution of the ability parameters. Default is 0 for each scale. +#' @param AbilitySD Numeric array. Standard deviation for the prior distribution of the ability parameters. Default is 10 for each scale. +#' @param AbilityCorr Matrix. Correlation matrix for the ability parameters. Default is an identity matrix. +#' @param AMeanSD Numeric. Standard deviation for the prior distribution of the discrimination parameters. Default is 1. +#' @param BMeanSD Numeric. Standard deviation for the prior distribution of the difficulty parameters. Default is \code{BSD}. +#' @param logitCMeanSD Numeric. Standard deviation for the prior distribution of the guessing parameters (on logit scale). Default is \code{logitCSD}. +#' @param logitDMeanSD Numeric. Standard deviation for the prior distribution of the upper asymptote parameters (on logit scale). Default is \code{logitDSD}. +#' @param AbilityMeanSD Numeric array. Standard deviation for the prior distribution of the ability parameters. Default is 1 for each scale. +#' @param iter Integer. Maximum number of iterations for the fitting algorithm. Default is 2000. +#' @param cores Integer. Number of cores to use for parallel computation. Default is 6. +#' @param carefulfit Logical. Whether to use a slower, careful fitting procedure. Default is FALSE. Experimental. +#' @param ebayes Logical. Whether to use empirical Bayes estimation. Default is TRUE. With ebayes, the priors are adapted based on a first pass estimate. +#' @param ebayesmultiplier Numeric. Multiplier for the widths of the empirical Bayes priors. Default is 2, as this appears to work better in practice. +#' @param ebayesFromFixed Logical. Whether to initialize empirical Bayes from any specifed fixed values Default is FALSE. +#' @param ebayesiter Integer. Number of iterations for empirical Bayes estimation. Default is 1. +#' @param estMeans Character vector. Which means to estimate from 'ability', 'A', 'B', 'C', 'D'. Default is c('ability', 'B', 'C', 'D'), with +#' discrimination means fixed. +#' @param priors Logical. Whether to use prior distributions. Default is TRUE. +#' @param integrateEachAbility Logical. Whether to integrate across each ability. Default is FALSE. +#' @param integrateEachAbilityFixedSE Logical. Whether to integrate each ability with fixed standard error. Default is FALSE. +#' @param mml Logical. Experimental and not working well. Whether to use marginal maximum likelihood estimation. Default is FALSE. +#' @param NintegratePoints Integer. Number of integration points for numerical integration. Default is 5. +#' @param normalise Logical. Whether to normalize the output estimates. Default is FALSE. +#' @param normaliseScale Numeric. Scale for normalization. Default is 1. +#' @param normaliseMean Numeric. Mean for normalization. Default is 0. +#' @param dropPerfectScores Logical. Whether to drop perfect scores from each subject and item before estimation. Default is TRUE. +#' @param trainingRows Integer vector. Rows of data to use for estimation of parameters. Default is all rows in \code{dat}. +#' @param init Initial values for the fitting algorithm. Default is NA. +#' @param tol Numeric. Tolerance for convergence. Default attempts to sensibly adjust for amount of data. +#' @param ... Additional arguments passed to the fitting function. +#' +#' @return A list containing the fitted IRT model parameters and additional information about the fit. #' -#' @param dat -#' @param score -#' @param id -#' @param item -#' @param scale -#' @param pl -#' @param Adata -#' @param Bdata -#' @param Cdata -#' @param personDat -#' @param invspAMeandat -#' @param invspASD -#' @param BMeandat -#' @param BSD -#' @param logitCMeandat -#' @param logitCSD -#' @param logitDMeandat -#' @param logitDSD -#' @param AbilityMeandat -#' @param AbilitySD -#' @param AMeanSD -#' @param BMeanSD -#' @param logitCMeanSD -#' @param logitDMeanSD -#' @param AbilityMeanSD -#' @param iter -#' @param cores -#' @param carefulfit -#' @param ebayes -#' @param ebayesmultiplier -#' @param estMeans -#' @param priors -#' @param normalise -#' @param trainingRows -#' @param tol -#' @param ... #' #' @return #' @export @@ -373,19 +396,21 @@ fitIRT <- function(dat,score='score', id='id', item='Item', scale='Scale',pl=1, AMeanSD=1,BMeanSD=BSD,logitCMeanSD=logitCSD,logitDMeanSD=logitDSD, AbilityMeanSD=array(1,dim=c(length(unique(dat[[scale]])))), iter=2000,cores=6,carefulfit=FALSE, - ebayes=TRUE,ebayesmultiplier=1,ebayesFromFixed=FALSE,ebayesiter=1, + ebayes=TRUE,ebayesmultiplier=2,ebayesFromFixed=FALSE,ebayesiter=1, estMeans=c('ability','B','C','D'),priors=TRUE, integrateEachAbility=FALSE, integrateEachAbilityFixedSE=FALSE, mml=FALSE,NintegratePoints=5, normalise=FALSE,normaliseScale=1,normaliseMean=0, dropPerfectScores=TRUE,trainingRows=1:nrow(dat), - init=NA,tol=1e-2,...){ + init=NA,tol=1e-8 * 10^(log(nrow(dat), 10)),...){ sdat <-list() #initialize standata object basetol=tol itemPreds <- unique(c(AitemPreds,BitemPreds,CitemPreds,DitemPreds)) + if(mml) stop('MML currently disabled -- hopefully working one day!') + #setup unlikely names to use in data.table calls to avoid overlap from user defined names idref. <- id; scaleref. <- scale; itemref. <- item; scoreref. <- score; @@ -723,7 +748,7 @@ fitIRT <- function(dat,score='score', id='id', item='Item', scale='Scale',pl=1, rownames(fit$pars$C)[itemIndex$new] <- itemIndex$original rownames(fit$pars$D)[itemIndex$new] <- itemIndex$original - return(fit) + return(list(fit=fit,sdat=sdat)) } rm(dat) @@ -741,6 +766,8 @@ fitIRT <- function(dat,score='score', id='id', item='Item', scale='Scale',pl=1, fit = fit, narrowPriors = JMLseq[[i]]$narrowPriors, tol=tol,...) + sdat <- fit$sdat + fit <- fit$fit } # if(ebayes) fit$fitML <- fitML diff --git a/bigIRT.Rproj b/bigIRT.Rproj index d40bd6e..6ae45e7 100644 --- a/bigIRT.Rproj +++ b/bigIRT.Rproj @@ -17,7 +17,6 @@ StripTrailingWhitespace: Yes BuildType: Package PackageUseDevtools: Yes -PackageCleanBeforeInstall: No PackageInstallArgs: --no-multiarch --with-keep.source PackageBuildArgs: --with-keep.source PackageRoxygenize: rd,collate,namespace diff --git a/inst/stan/irtmml.stan b/inst/stan/irtmml.bak similarity index 100% rename from inst/stan/irtmml.stan rename to inst/stan/irtmml.bak diff --git a/readme.Rmd b/readme.Rmd index 7c532a3..4e49e66 100644 --- a/readme.Rmd +++ b/readme.Rmd @@ -1,7 +1,7 @@ --- title: "bigIRT" author: "Charles Driver" -date: "12/05/2021" +date: "" output: github_document --- @@ -13,6 +13,9 @@ knitr::opts_chunk$set(echo = TRUE) ``` ## Install +If building on windows, first ensure you have Rtools installed. https://cran.r-project.org/bin/windows/Rtools/ + +Run the following code to install the package: ```{r cars, eval=FALSE} remotes::install_github('cdriveraus/bigIRT', INSTALL_opts = "--no-multiarch", dependencies = c("Depends", "Imports")) @@ -20,33 +23,80 @@ remotes::install_github('cdriveraus/bigIRT', INSTALL_opts = "--no-multiarch", de ## Example +This example generates data from a 3 parameter logistic model, fits the model using mirt and bigIRT, and compares the parameter estimates. + ```{r, results=FALSE,plot=FALSE} library(bigIRT) -#Generate some data (here 2pl model -require(data.table) -dat <- simIRT(Nsubs = 1000,Nitems = 100,Nscales = 1, - logitCMean = -10,logitCSD = 0,AMean = 1,ASD = .3, +library(data.table) +library(mirt) +library(ggplot2) + +#generate data +dat <- simIRT(Nsubs = 1000,Nitems = 50,Nscales = 1, + logitCMean = -2,logitCSD = 1,AMean = 1,ASD = .3, BMean=0,BSD = .5, AbilityMean = 0,AbilitySD = 1) -#convert to wide for TAM +#normalise true parameters and store +trueitempars <- normaliseIRT(B = dat$B, A = dat$A, Ability = dat$Ability) + +#convert to wide for mirt wdat <- data.frame(dcast(data.table(dat$dat),formula = 'id ~ Item',value.var='score')[,-1]) -#fit using TAM -require(TAM) -tfit <-tam.mml.2pl(resp = wdat,est.variance = TRUE) +#fit using mirt +mirtfit <- mirt(wdat,model=1, itemtype = '3PL',method = "EM") +#get normalized parameter estimates +mirtItempars <- coef(mirtfit) +mirtItempars$GroupPars <- NULL +mirtItempars <- rbindlist(lapply(mirtItempars, function(x) as.data.frame(x))) +mirtPersonpars <- fscores(mirtfit)[,1] +mirtpars<- normaliseIRT(B = -mirtItempars$d / mirtItempars$a1, A = mirtItempars$a1, Ability = mirtPersonpars) #fit using bigIRT -fit <- fitIRT(dat$dat,cores=2,score = 'score',id = 'id',scale = 'Scale',item = 'Item',pl=2,verbose=0,plot=0) +bigirtfit <- fitIRT(dat$dat, pl=3, cores=4, dropPerfectScores = FALSE, stochastic=TRUE, + score = 'score', id = 'id', scale = 'Scale', item = 'Item') + +#get normalized parameter estimates +bigirtpars <- normaliseIRT(B = bigirtfit$itemPars$B, A = bigirtfit$itemPars$A, Ability = bigirtfit$pars$Ability) + +#combine parameter outputs +itempars <- rbind( + data.table(Estimator = 'mirt', Parameter = 'A', Estimate=c(mirtpars$A), TrueValue=c(trueitempars$A)), + data.table(Estimator = 'mirt', Parameter = 'B', Estimate=c(mirtpars$B), TrueValue=c(trueitempars$B)), + data.table(Estimator = 'mirt', Parameter = 'C', Estimate=c(mirtItempars$g), TrueValue=c(dat$C)), + data.table(Estimator = 'bigIRT', Parameter = 'A', Estimate=c(bigirtpars$A), TrueValue=c(trueitempars$A)), + data.table(Estimator = 'bigIRT', Parameter = 'B', Estimate=c(bigirtpars$B), TrueValue=c(trueitempars$B)), + data.table(Estimator = 'bigIRT', Parameter = 'C', Estimate=c(bigirtfit$itemPars$C), TrueValue=c(dat$C)) +) + +personpars <- rbind( + data.table(Estimator = 'mirt', Parameter = 'Ability', Estimate=c(mirtpars$Ability), TrueValue=c(dat$Ability)), + data.table(Estimator = 'bigIRT', Parameter = 'Ability', Estimate=c(bigirtpars$Ability), TrueValue=c(dat$Ability)) +) + +pars <- rbind(itempars,personpars) + +#plot parameter estimates +ggplot(pars,aes(x=TrueValue,y=Estimate-TrueValue,color=Estimator)) + + geom_point() + + facet_wrap(~Parameter,scales='free')+ + theme_bw()+ + geom_hline(yintercept=0,linetype='dashed') + +#plot average absolute error +ggplot(pars,aes(x=TrueValue,y=abs(Estimate-TrueValue),colour=Estimator, fill=Estimator))+ + geom_smooth()+ + facet_wrap(~Parameter,scales='free')+ + theme_bw()+ + coord_cartesian(ylim=c(0,NA))+ + geom_hline(yintercept=0,linetype='dashed') -#some summary stuff: -plot(dat$Ability,(fit$pars$Ability-dat$Ability)^2) #Ability error given Ability -sqrt(mean((fit$pars$Ability-dat$Ability)^2)) #rms error stat +#plot average error +ggplot(pars,aes(x=TrueValue,y=Estimate-TrueValue,colour=Estimator, fill=Estimator))+ + geom_smooth()+ + facet_wrap(~Parameter,scales='free')+ + theme_bw()+geom_hline(yintercept=0,linetype='dashed') -#correlations of estimated vs true -cor(data.frame(True=dat$Ability,Est=fit$pars$Ability)) -cor(data.frame(True=dat$A,Est=fit$pars$A)) -cor(data.frame(True=dat$B,Est=fit$pars$B)) ``` diff --git a/readme.md b/readme.md index 0707c8d..8d1d6b4 100644 --- a/readme.md +++ b/readme.md @@ -1,7 +1,6 @@ bigIRT ================ Charles Driver -12/05/2021 [![R-CMD-check](https://github.com/cdriveraus/bigIRT/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/cdriveraus/bigIRT/actions/workflows/R-CMD-check.yaml) @@ -9,73 +8,131 @@ Charles Driver ## Install +If building on windows, first ensure you have Rtools installed. + + +Run the following code to install the package: + ``` r remotes::install_github('cdriveraus/bigIRT', INSTALL_opts = "--no-multiarch", dependencies = c("Depends", "Imports")) ``` ## Example +This example generates data from a 3 parameter logistic model, fits the +model using mirt and bigIRT, and compares the parameter estimates. + ``` r library(bigIRT) ``` ## Loading required package: data.table - ## Warning: package 'data.table' was built under R version 4.2.2 +``` r +library(data.table) +library(mirt) +``` + + ## Loading required package: stats4 + + ## Loading required package: lattice ``` r -#Generate some data (here 2pl model -require(data.table) -dat <- simIRT(Nsubs = 1000,Nitems = 100,Nscales = 1, - logitCMean = -10,logitCSD = 0,AMean = 1,ASD = .3, +library(ggplot2) + +#generate data +dat <- simIRT(Nsubs = 1000,Nitems = 50,Nscales = 1, + logitCMean = -2,logitCSD = 1,AMean = 1,ASD = .3, BMean=0,BSD = .5, AbilityMean = 0,AbilitySD = 1) -#convert to wide for TAM +#normalise true parameters and store +trueitempars <- normaliseIRT(B = dat$B, A = dat$A, Ability = dat$Ability) + +#convert to wide for mirt wdat <- data.frame(dcast(data.table(dat$dat),formula = 'id ~ Item',value.var='score')[,-1]) -#fit using TAM -require(TAM) +#fit using mirt +mirtfit <- mirt(wdat,model=1, itemtype = '3PL',method = "EM") ``` - ## Loading required package: TAM - - ## Loading required package: CDM - - ## Loading required package: mvtnorm - - ## ********************************** - ## ** CDM 8.2-6 (2022-08-25 15:43:23) - ## ** Cognitive Diagnostic Models ** - ## ********************************** - - ## * TAM 4.1-4 (2022-08-28 16:03:54) + ## EM cycles terminated after 500 iterations. ``` r -tfit <-tam.mml.2pl(resp = wdat,est.variance = TRUE) - +#get normalized parameter estimates +mirtItempars <- coef(mirtfit) +mirtItempars$GroupPars <- NULL +mirtItempars <- rbindlist(lapply(mirtItempars, function(x) as.data.frame(x))) +mirtPersonpars <- fscores(mirtfit)[,1] +mirtpars<- normaliseIRT(B = -mirtItempars$d / mirtItempars$a1, A = mirtItempars$a1, Ability = mirtPersonpars) #fit using bigIRT -fit <- fitIRT(dat$dat,cores=2,score = 'score',id = 'id',scale = 'Scale',item = 'Item',pl=2,verbose=0,plot=0) +bigirtfit <- fitIRT(dat$dat, pl=3, cores=4, dropPerfectScores = FALSE, stochastic=TRUE, + score = 'score', id = 'id', scale = 'Scale', item = 'Item') ``` ## Free estimation step... + ## Progress est. = 0%, LPchange = 96, Iter = 31, LP = -31881.183848043Progress est. = 0%, LPchange = 120, Iter = 32, LP = -30989.9833787027Progress est. = 0%, LPchange = 140, Iter = 33, LP = -30315.1728567025Progress est. = 0.15%, LPchange = 140, Iter = 34, LP = -30315.1728567025Progress est. = 0%, LPchange = 160, Iter = 35, LP = -29758.8949673376Progress est. = 0.07%, LPchange = 160, Iter = 36, LP = -29741.7023523809Progress est. = 0.17%, LPchange = 150, Iter = 37, LP = -29741.7023523809Progress est. = 0%, LPchange = 160, Iter = 38, LP = -29540.1019415843Progress est. = 0.08%, LPchange = 160, Iter = 39, LP = -29540.1019415843Progress est. = 0%, LPchange = 160, Iter = 40, LP = -29364.731277841Progress est. = 0.08%, LPchange = 160, Iter = 41, LP = -29364.731277841Progress est. = 0%, LPchange = 160, Iter = 42, LP = -29282.2030974931Progress est. = 0.05%, LPchange = 160, Iter = 43, LP = -29271.1779731264Progress est. = 0.12%, LPchange = 160, Iter = 44, LP = -29271.1779731264Progress est. = 0%, LPchange = 160, Iter = 45, LP = -29167.404110647Progress est. = 0.07%, LPchange = 160, Iter = 46, LP = -29167.404110647Progress est. = 0.05%, LPchange = 160, Iter = 47, LP = -29125.2752185988Progress est. = 0.12%, LPchange = 160, Iter = 48, LP = -29125.2752185988Progress est. = 0.07%, LPchange = 160, Iter = 49, LP = -29067.189532407Progress est. = 0.14%, LPchange = 160, Iter = 50, LP = -29067.189532407Progress est. = 0.15%, LPchange = 160, Iter = 51, LP = -29042.7679410761Progress est. = 0.25%, LPchange = 160, Iter = 52, LP = -29042.7679410761Progress est. = 0.35%, LPchange = 160, Iter = 53, LP = -29019.2466671163Progress est. = 0.56%, LPchange = 150, Iter = 54, LP = -29019.2466671163Progress est. = 0.86%, LPchange = 150, Iter = 55, LP = -29005.8311343831Progress est. = 1.33%, LPchange = 140, Iter = 56, LP = -29005.8311343831Progress est. = 1.98%, LPchange = 130, Iter = 57, LP = -28992.514047427Progress est. = 2.88%, LPchange = 120, Iter = 58, LP = -28992.514047427Progress est. = 3.86%, LPchange = 110, Iter = 59, LP = -28982.45257715Progress est. = 4.54%, LPchange = 100, Iter = 60, LP = -28982.45257715Progress est. = 5.01%, LPchange = 97, Iter = 61, LP = -28975.9851289883Progress est. = 8.53%, LPchange = 67, Iter = 62, LP = -28975.9851289883Progress est. = 12.42%, LPchange = 45, Iter = 63, LP = -28970.4040254059Progress est. = 12.42%, LPchange = 45, Iter = 64, LP = -28970.4040254059Progress est. = 17.55%, LPchange = 26, Iter = 65, LP = -28969.8312804693Progress est. = 17.76%, LPchange = 26, Iter = 66, LP = -28969.8312804693Progress est. = 17.74%, LPchange = 26, Iter = 67, LP = -28967.8026396441Progress est. = 20.64%, LPchange = 19, Iter = 68, LP = -28967.8026396441Progress est. = 20.58%, LPchange = 19, Iter = 69, LP = -28964.2404741661Progress est. = 24.08%, LPchange = 13, Iter = 70, LP = -28964.2404741661Progress est. = 24.01%, LPchange = 13, Iter = 71, LP = -28961.4045710163Progress est. = 26.21%, LPchange = 11, Iter = 72, LP = -28961.4045710163Progress est. = 26.53%, LPchange = 10, Iter = 73, LP = -28960.8663423145Progress est. = 26.53%, LPchange = 10, Iter = 74, LP = -28960.8663423145Progress est. = 30.27%, LPchange = 7, Iter = 75, LP = -28957.0685646598Progress est. = 30.27%, LPchange = 7, Iter = 76, LP = -28957.0685646598Progress est. = 32.31%, LPchange = 5.7, Iter = 77, LP = -28955.067960893Progress est. = 32.31%, LPchange = 5.7, Iter = 78, LP = -28955.067960893Progress est. = 36.26%, LPchange = 3.8, Iter = 79, LP = -28954.2453771907Progress est. = 36.26%, LPchange = 3.8, Iter = 80, LP = -28954.2453771907Progress est. = 38.49%, LPchange = 3, Iter = 81, LP = -28953.1914520241Progress est. = 38.49%, LPchange = 3, Iter = 82, LP = -28953.1914520241Progress est. = 41.35%, LPchange = 2.2, Iter = 83, LP = -28952.6949215916Progress est. = 41.35%, LPchange = 2.2, Iter = 84, LP = -28952.6949215916Progress est. = 43.4%, LPchange = 1.8, Iter = 85, LP = -28952.0919978326Progress est. = 43.4%, LPchange = 1.8, Iter = 86, LP = -28952.0919978326Progress est. = 45.86%, LPchange = 1.4, Iter = 87, LP = -28950.8835694868Progress est. = 45.86%, LPchange = 1.4, Iter = 88, LP = -28950.8835694868Progress est. = 48.44%, LPchange = 1.1, Iter = 89, LP = -28950.6164306822Progress est. = 48.44%, LPchange = 1.1, Iter = 90, LP = -28950.6164306822Progress est. = 50.13%, LPchange = 0.89, Iter = 91, LP = -28949.2621311744Progress est. = 50.13%, LPchange = 0.89, Iter = 92, LP = -28949.2621311744Progress est. = 52.29%, LPchange = 0.71, Iter = 93, LP = -28949.0684669031Progress est. = 52.29%, LPchange = 0.71, Iter = 94, LP = -28949.0684669031Progress est. = 52.27%, LPchange = 0.71, Iter = 95, LP = -28948.4442742876Progress est. = 52.27%, LPchange = 0.71, Iter = 96, LP = -28948.4442742876Progress est. = 53.04%, LPchange = 0.66, Iter = 97, LP = -28948.0518640552Progress est. = 53.04%, LPchange = 0.66, Iter = 98, LP = -28948.0518640552Progress est. = 54.73%, LPchange = 0.55, Iter = 99, LP = -28947.685346405Progress est. = 54.73%, LPchange = 0.55, Iter = 100, LP = -28947.685346405Progress est. = 56.54%, LPchange = 0.46, Iter = 101, LP = -28947.685346405Progress est. = 56.39%, LPchange = 0.46, Iter = 102, LP = -28947.4699300992Progress est. = 56.73%, LPchange = 0.45, Iter = 103, LP = -28947.4109751839Progress est. = 56.73%, LPchange = 0.45, Iter = 104, LP = -28947.4109751839Progress est. = 59.37%, LPchange = 0.34, Iter = 105, LP = -28946.8404582483Progress est. = 59.33%, LPchange = 0.34, Iter = 106, LP = -28946.8039024152Progress est. = 61.29%, LPchange = 0.28, Iter = 107, LP = -28946.6906223235Progress est. = 61.29%, LPchange = 0.28, Iter = 108, LP = -28946.6906223235Progress est. = 61.82%, LPchange = 0.26, Iter = 109, LP = -28946.3188242225Progress est. = 61.8%, LPchange = 0.26, Iter = 110, LP = -28946.3027848838Progress est. = 62.63%, LPchange = 0.24, Iter = 111, LP = -28945.9000603619Progress est. = 62.63%, LPchange = 0.24, Iter = 112, LP = -28945.9000603619Progress est. = 63.3%, LPchange = 0.23, Iter = 113, LP = -28945.9000603619Progress est. = 63.19%, LPchange = 0.23, Iter = 114, LP = -28945.8207003791Progress est. = 64.08%, LPchange = 0.21, Iter = 115, LP = -28945.8207003791Progress est. = 63.4%, LPchange = 0.22, Iter = 116, LP = -28945.3618547292Progress est. = 65.3%, LPchange = 0.18, Iter = 117, LP = -28945.3618547292Progress est. = 64.03%, LPchange = 0.21, Iter = 118, LP = -28944.5789290077Progress est. = 64.44%, LPchange = 0.2, Iter = 119, LP = -28944.5789290077Progress est. = 63.3%, LPchange = 0.23, Iter = 120, LP = -28943.8192168316Progress est. = 65.44%, LPchange = 0.18, Iter = 121, LP = -28943.8192168316Progress est. = 64.3%, LPchange = 0.2, Iter = 122, LP = -28943.1354249754Progress est. = 64.61%, LPchange = 0.2, Iter = 123, LP = -28943.1354249754Progress est. = 62.52%, LPchange = 0.25, Iter = 124, LP = -28941.6958772409Progress est. = 63.37%, LPchange = 0.22, Iter = 125, LP = -28941.6958772409Progress est. = 61.02%, LPchange = 0.29, Iter = 126, LP = -28939.8283147545Progress est. = 61.47%, LPchange = 0.27, Iter = 127, LP = -28939.8283147545Progress est. = 59.75%, LPchange = 0.33, Iter = 128, LP = -28938.2257491559Progress est. = 60.12%, LPchange = 0.32, Iter = 129, LP = -28938.2257491559Progress est. = 59.43%, LPchange = 0.34, Iter = 130, LP = -28937.5268518013Progress est. = 58.72%, LPchange = 0.36, Iter = 131, LP = -28936.7485062605Progress est. = 58.92%, LPchange = 0.36, Iter = 132, LP = -28936.7485062605Progress est. = 58.28%, LPchange = 0.38, Iter = 133, LP = -28935.9516154424Progress est. = 58.16%, LPchange = 0.39, Iter = 134, LP = -28935.8100837686Progress est. = 58.17%, LPchange = 0.39, Iter = 135, LP = -28935.2595161092Progress est. = 57.83%, LPchange = 0.4, Iter = 136, LP = -28934.8020477989Progress est. = 57.75%, LPchange = 0.4, Iter = 137, LP = -28934.5849597471Progress est. = 57.52%, LPchange = 0.41, Iter = 138, LP = -28934.2899134252Progress est. = 57.4%, LPchange = 0.42, Iter = 139, LP = -28933.7724121584Progress est. = 57.32%, LPchange = 0.42, Iter = 140, LP = -28933.6530269976Progress est. = 57.36%, LPchange = 0.42, Iter = 141, LP = -28933.3020278109Progress est. = 56.98%, LPchange = 0.44, Iter = 142, LP = -28932.795121071Progress est. = 56.89%, LPchange = 0.44, Iter = 143, LP = -28932.664859981Progress est. = 56.6%, LPchange = 0.45, Iter = 144, LP = -28932.181230873Progress est. = 56.44%, LPchange = 0.46, Iter = 145, LP = -28931.9535752563Progress est. = 56.46%, LPchange = 0.46, Iter = 146, LP = -28931.5203962935Progress est. = 56.17%, LPchange = 0.48, Iter = 147, LP = -28931.1014560262Progress est. = 56.67%, LPchange = 0.45, Iter = 148, LP = -28931.0336445755Progress est. = 56.51%, LPchange = 0.46, Iter = 149, LP = -28930.8087410397Progress est. = 57%, LPchange = 0.44, Iter = 150, LP = -28930.7340342276Progress est. = 56.37%, LPchange = 0.47, Iter = 151, LP = -28929.8485621478Progress est. = 56.32%, LPchange = 0.47, Iter = 152, LP = -28929.1000140303Progress est. = 56.32%, LPchange = 0.47, Iter = 153, LP = -28929.1000140303Progress est. = 57.01%, LPchange = 0.44, Iter = 154, LP = -28928.6238941825Progress est. = 57.01%, LPchange = 0.44, Iter = 155, LP = -28928.6238941825Progress est. = 57.93%, LPchange = 0.4, Iter = 156, LP = -28927.9503713177Progress est. = 57.93%, LPchange = 0.4, Iter = 157, LP = -28927.9503713177Progress est. = 59.25%, LPchange = 0.35, Iter = 158, LP = -28927.8664724846Progress est. = 58.66%, LPchange = 0.37, Iter = 159, LP = -28927.2119747574Progress est. = 58.84%, LPchange = 0.36, Iter = 160, LP = -28926.7188699041Progress est. = 58.77%, LPchange = 0.36, Iter = 161, LP = -28925.8628408302Progress est. = 58.4%, LPchange = 0.38, Iter = 162, LP = -28925.4387203277Progress est. = 58.71%, LPchange = 0.37, Iter = 163, LP = -28924.9988512583Progress est. = 58.38%, LPchange = 0.38, Iter = 164, LP = -28924.4704831714Progress est. = 58.31%, LPchange = 0.38, Iter = 165, LP = -28923.8362586693Progress est. = 58.52%, LPchange = 0.37, Iter = 166, LP = -28923.633297779Progress est. = 58.31%, LPchange = 0.38, Iter = 167, LP = -28923.1718189851Progress est. = 58.23%, LPchange = 0.38, Iter = 168, LP = -28922.7755633244Progress est. = 58.38%, LPchange = 0.38, Iter = 169, LP = -28922.4367577876Progress est. = 58.24%, LPchange = 0.38, Iter = 170, LP = -28922.1536817172Progress est. = 58.26%, LPchange = 0.38, Iter = 171, LP = -28921.8223411671Progress est. = 58.48%, LPchange = 0.37, Iter = 172, LP = -28921.5743037074Progress est. = 58.36%, LPchange = 0.38, Iter = 173, LP = -28921.3103253033Progress est. = 58.52%, LPchange = 0.37, Iter = 174, LP = -28921.0073008899Progress est. = 58.5%, LPchange = 0.37, Iter = 175, LP = -28920.7598231573Progress est. = 58.69%, LPchange = 0.37, Iter = 176, LP = -28920.5469887493Progress est. = 58.81%, LPchange = 0.36, Iter = 177, LP = -28920.2657072933Progress est. = 58.68%, LPchange = 0.37, Iter = 178, LP = -28920.0469825412Progress est. = 58.74%, LPchange = 0.36, Iter = 179, LP = -28919.8895480458Progress est. = 58.81%, LPchange = 0.36, Iter = 180, LP = -28919.8895480458Progress est. = 59.63%, LPchange = 0.33, Iter = 181, LP = -28919.8895480458Progress est. = 60.38%, LPchange = 0.31, Iter = 182, LP = -28919.8895480458Progress est. = 60.17%, LPchange = 0.31, Iter = 183, LP = -28919.6883852369Progress est. = 60.16%, LPchange = 0.31, Iter = 184, LP = -28919.2028816212Progress est. = 60.16%, LPchange = 0.31, Iter = 185, LP = -28919.2028816212Progress est. = 60.87%, LPchange = 0.29, Iter = 186, LP = -28919.2028816212Progress est. = 60.41%, LPchange = 0.31, Iter = 187, LP = -28918.7741768181Progress est. = 60.46%, LPchange = 0.3, Iter = 188, LP = -28918.7387203504Progress est. = 60.86%, LPchange = 0.29, Iter = 189, LP = -28918.4505081283Progress est. = 61.27%, LPchange = 0.28, Iter = 190, LP = -28918.3236001957Progress est. = 62.03%, LPchange = 0.26, Iter = 191, LP = -28918.1058073992Progress est. = 62.51%, LPchange = 0.25, Iter = 192, LP = -28918.0570317659Progress est. = 62.89%, LPchange = 0.24, Iter = 193, LP = -28917.9049684733Progress est. = 63.49%, LPchange = 0.22, Iter = 194, LP = -28917.8037214146Progress est. = 64.35%, LPchange = 0.2, Iter = 195, LP = -28917.7406046653Progress est. = 64.52%, LPchange = 0.2, Iter = 196, LP = -28917.6434591612Progress est. = 65.1%, LPchange = 0.19, Iter = 197, LP = -28917.5350478169Progress est. = 65.67%, LPchange = 0.18, Iter = 198, LP = -28917.4608338844Progress est. = 66.2%, LPchange = 0.17, Iter = 199, LP = -28917.4059806943Progress est. = 66.55%, LPchange = 0.16, Iter = 200, LP = -28917.3029514232Progress est. = 67.15%, LPchange = 0.15, Iter = 201, LP = -28917.2654708159Progress est. = 67.54%, LPchange = 0.15, Iter = 202, LP = -28917.2004127303Progress est. = 67.96%, LPchange = 0.14, Iter = 203, LP = -28917.1212326491Progress est. = 68.56%, LPchange = 0.13, Iter = 204, LP = -28917.0700207521Progress est. = 69.14%, LPchange = 0.12, Iter = 205, LP = -28917.0529493555Progress est. = 69.47%, LPchange = 0.12, Iter = 206, LP = -28916.9672329827Progress est. = 70.08%, LPchange = 0.11, Iter = 207, LP = -28916.9063531024Progress est. = 70.6%, LPchange = 0.11, Iter = 208, LP = -28916.8639995257Progress est. = 70.84%, LPchange = 0.1, Iter = 209, LP = -28916.7844885762Progress est. = 70.64%, LPchange = 0.11, Iter = 210, LP = -28916.7183173174Progress est. = 70.45%, LPchange = 0.11, Iter = 211, LP = -28916.6571189945Progress est. = 70.32%, LPchange = 0.11, Iter = 212, LP = -28916.6134571183Progress est. = 70.79%, LPchange = 0.1, Iter = 213, LP = -28916.5653876638Progress est. = 72.18%, LPchange = 0.09, Iter = 214, LP = -28916.5021541044Progress est. = 72.03%, LPchange = 0.092, Iter = 215, LP = -28916.4576663505Progress est. = 71.85%, LPchange = 0.093, Iter = 216, LP = -28916.4077440467Progress est. = 73.28%, LPchange = 0.08, Iter = 217, LP = -28916.3632722288Progress est. = 73.19%, LPchange = 0.081, Iter = 218, LP = -28916.3056911788Progress est. = 74.4%, LPchange = 0.071, Iter = 219, LP = -28916.3056911788Progress est. = 74.93%, LPchange = 0.068, Iter = 220, LP = -28916.2931030741Progress est. = 75.44%, LPchange = 0.064, Iter = 221, LP = -28916.1801901689Progress est. = 75.38%, LPchange = 0.065, Iter = 222, LP = -28916.1196856195Progress est. = 75.98%, LPchange = 0.061, Iter = 223, LP = -28916.0838567065Progress est. = 76.35%, LPchange = 0.058, Iter = 224, LP = -28916.0527835473Progress est. = 76.24%, LPchange = 0.059, Iter = 225, LP = -28915.969318866Progress est. = 76.59%, LPchange = 0.057, Iter = 226, LP = -28915.9347847458Progress est. = 77.22%, LPchange = 0.053, Iter = 227, LP = -28915.9347847458Progress est. = 76.97%, LPchange = 0.055, Iter = 228, LP = -28915.81792386Progress est. = 77.19%, LPchange = 0.054, Iter = 229, LP = -28915.8008011419Progress est. = 77.29%, LPchange = 0.053, Iter = 230, LP = -28915.7151493511Progress est. = 77.52%, LPchange = 0.052, Iter = 231, LP = -28915.7151493511Progress est. = 77.77%, LPchange = 0.05, Iter = 232, LP = -28915.6895034644Progress est. = 78.1%, LPchange = 0.049, Iter = 233, LP = -28915.6609702862Progress est. = 78.2%, LPchange = 0.048, Iter = 234, LP = -28915.6243054875Progress est. = 78.13%, LPchange = 0.049, Iter = 235, LP = -28915.5967315976Progress est. = 78.35%, LPchange = 0.047, Iter = 236, LP = -28915.5438662636Progress est. = 78.6%, LPchange = 0.046, Iter = 237, LP = -28915.5199515292Progress est. = 78.66%, LPchange = 0.046, Iter = 238, LP = -28915.4857740791Progress est. = 79.04%, LPchange = 0.044, Iter = 239, LP = -28915.4594247337Progress est. = 79.38%, LPchange = 0.043, Iter = 240, LP = -28915.4404355601Progress est. = 79.63%, LPchange = 0.042, Iter = 241, LP = -28915.4113466016Progress est. = 79.84%, LPchange = 0.041, Iter = 242, LP = -28915.3950055179Progress est. = 80.01%, LPchange = 0.04, Iter = 243, LP = -28915.3674987159Progress est. = 80.4%, LPchange = 0.038, Iter = 244, LP = -28915.3520246885Progress est. = 80.55%, LPchange = 0.038, Iter = 245, LP = -28915.3261270956Progress est. = 80.86%, LPchange = 0.037, Iter = 246, LP = -28915.3115795559Progress est. = 81.05%, LPchange = 0.036, Iter = 247, LP = -28915.2885553468Progress est. = 81.4%, LPchange = 0.035, Iter = 248, LP = -28915.2694477245Progress est. = 81.4%, LPchange = 0.035, Iter = 249, LP = -28915.2694477245Progress est. = 81.52%, LPchange = 0.034, Iter = 250, LP = -28915.2694477245Progress est. = 82.43%, LPchange = 0.031, Iter = 251, LP = -28915.2494377984Progress est. = 83.08%, LPchange = 0.029, Iter = 252, LP = -28915.2494377984Progress est. = 83.43%, LPchange = 0.028, Iter = 253, LP = -28915.2445221812Progress est. = 83.5%, LPchange = 0.028, Iter = 254, LP = -28915.2194589343Progress est. = 84.42%, LPchange = 0.025, Iter = 255, LP = -28915.2122095364Progress est. = 84.52%, LPchange = 0.025, Iter = 256, LP = -28915.1850214512Progress est. = 84.33%, LPchange = 0.025, Iter = 257, LP = -28915.170277358Progress est. = 85.82%, LPchange = 0.022, Iter = 258, LP = -28915.1630059837Progress est. = 85.83%, LPchange = 0.022, Iter = 259, LP = -28915.1466377129Progress est. = 87.04%, LPchange = 0.019, Iter = 260, LP = -28915.1381780003Progress est. = 86.81%, LPchange = 0.02, Iter = 261, LP = -28915.1247022068Progress est. = 87.1%, LPchange = 0.019, Iter = 262, LP = -28915.1166256806Progress est. = 87.47%, LPchange = 0.018, Iter = 263, LP = -28915.1096820283Progress est. = 87.85%, LPchange = 0.018, Iter = 264, LP = -28915.0942644383Progress est. = 88.2%, LPchange = 0.017, Iter = 265, LP = -28915.0852445227Progress est. = 89.09%, LPchange = 0.016, Iter = 266, LP = -28915.0775979615Progress est. = 89.55%, LPchange = 0.015, Iter = 267, LP = -28915.0755865146Progress est. = 90.06%, LPchange = 0.014, Iter = 268, LP = -28915.0641984226Progress est. = 90.44%, LPchange = 0.013, Iter = 269, LP = -28915.05450836Progress est. = 90.8%, LPchange = 0.013, Iter = 270, LP = -28915.0501702976Progress est. = 91.54%, LPchange = 0.012, Iter = 271, LP = -28915.0501702976Progress est. = 91.99%, LPchange = 0.011, Iter = 272, LP = -28915.0501702976Progress est. = 92.79%, LPchange = 0.011, Iter = 273, LP = -28915.0501702976Progress est. = 92.91%, LPchange = 0.01, Iter = 274, LP = -28915.0387920581Progress est. = 93.68%, LPchange = 0.0096, Iter = 275, LP = -28915.0369256734Progress est. = 93.71%, LPchange = 0.0096, Iter = 276, LP = -28915.0231089041Progress est. = 94.25%, LPchange = 0.0091, Iter = 277, LP = -28915.0159254562Progress est. = 94.79%, LPchange = 0.0086, Iter = 278, LP = -28915.0116616446Progress est. = 94.56%, LPchange = 0.0088, Iter = 279, LP = -28915.0054430346Progress est. = 94.38%, LPchange = 0.009, Iter = 280, LP = -28915.0004861775Progress est. = 94.96%, LPchange = 0.0084, Iter = 281, LP = -28914.996214417Progress est. = 94.83%, LPchange = 0.0086, Iter = 282, LP = -28914.9927473509Progress est. = 94.89%, LPchange = 0.0085, Iter = 283, LP = -28914.9895228775Progress est. = 95.74%, LPchange = 0.0078, Iter = 284, LP = -28914.986040929Progress est. = 95.93%, LPchange = 0.0076, Iter = 285, LP = -28914.9832666893Progress est. = 97.06%, LPchange = 0.0068, Iter = 286, LP = -28914.9814370904Progress est. = 97.66%, LPchange = 0.0064, Iter = 287, LP = -28914.9790448186Progress est. = 97.97%, LPchange = 0.0062, Iter = 288, LP = -28914.9777721629Progress est. = 98.86%, LPchange = 0.0056, Iter = 289, LP = -28914.9777721629Progress est. = 99.35%, LPchange = 0.0053, Iter = 290, LP = -28914.9777721629Progress est. = 100%, LPchange = 0.0049, Iter = 291, LP = -28914.9777721629 + ## Converged -- lp change within tol(0.005) ## Empirical Bayes step... + ## Progress est. = 0%, LPchange = 86, Iter = 31, LP = -32123.826767118Progress est. = 0%, LPchange = 110, Iter = 32, LP = -31332.3812540841Progress est. = 0%, LPchange = 140, Iter = 33, LP = -30331.4805795993Progress est. = 0.13%, LPchange = 140, Iter = 34, LP = -30331.4805795993Progress est. = 0%, LPchange = 140, Iter = 35, LP = -30177.6408870137Progress est. = 0.1%, LPchange = 140, Iter = 36, LP = -30177.6408870137Progress est. = 0%, LPchange = 140, Iter = 37, LP = -29965.8417245009Progress est. = 0.08%, LPchange = 140, Iter = 38, LP = -29965.8417245009Progress est. = 0.02%, LPchange = 140, Iter = 39, LP = -29890.9174012707Progress est. = 0.1%, LPchange = 140, Iter = 40, LP = -29890.9174012707Progress est. = 0%, LPchange = 150, Iter = 41, LP = -29684.7901678799Progress est. = 0.07%, LPchange = 150, Iter = 42, LP = -29684.7901678799Progress est. = 0%, LPchange = 150, Iter = 43, LP = -29548.4621545618Progress est. = 0.06%, LPchange = 150, Iter = 44, LP = -29548.4621545618Progress est. = 0%, LPchange = 150, Iter = 45, LP = -29467.1241605572Progress est. = 0.06%, LPchange = 150, Iter = 46, LP = -29467.1241605572Progress est. = 0.02%, LPchange = 150, Iter = 47, LP = -29406.2118255273Progress est. = 0.08%, LPchange = 150, Iter = 48, LP = -29406.2118255273Progress est. = 0.07%, LPchange = 150, Iter = 49, LP = -29370.7511833959Progress est. = 0.13%, LPchange = 150, Iter = 50, LP = -29370.7511833959Progress est. = 0.11%, LPchange = 150, Iter = 51, LP = -29325.2426426557Progress est. = 0.2%, LPchange = 150, Iter = 52, LP = -29325.2426426557Progress est. = 0.25%, LPchange = 150, Iter = 53, LP = -29282.5212426137Progress est. = 0.45%, LPchange = 140, Iter = 54, LP = -29282.5212426137Progress est. = 0.7%, LPchange = 140, Iter = 55, LP = -29259.5288526293Progress est. = 1.13%, LPchange = 130, Iter = 56, LP = -29259.5288526293Progress est. = 1.66%, LPchange = 120, Iter = 57, LP = -29221.2962544466Progress est. = 2.45%, LPchange = 110, Iter = 58, LP = -29221.2962544466Progress est. = 3.17%, LPchange = 100, Iter = 59, LP = -29201.8853365323Progress est. = 3.41%, LPchange = 98, Iter = 60, LP = -29201.8853365323Progress est. = 3.4%, LPchange = 98, Iter = 61, LP = -29184.8649625834Progress est. = 5.88%, LPchange = 72, Iter = 62, LP = -29184.8649625834Progress est. = 10.72%, LPchange = 39, Iter = 63, LP = -29164.3317549016Progress est. = 10.7%, LPchange = 39, Iter = 64, LP = -29161.6333097946Progress est. = 11.77%, LPchange = 34, Iter = 65, LP = -29156.0700863159Progress est. = 11.68%, LPchange = 34, Iter = 66, LP = -29143.5756485601Progress est. = 13.4%, LPchange = 28, Iter = 67, LP = -29133.758492946Progress est. = 13.36%, LPchange = 28, Iter = 68, LP = -29129.5737940532Progress est. = 13.99%, LPchange = 26, Iter = 69, LP = -29118.4315085524Progress est. = 13.97%, LPchange = 26, Iter = 70, LP = -29116.8368704746Progress est. = 16.31%, LPchange = 19, Iter = 71, LP = -29108.5459496927Progress est. = 16.18%, LPchange = 20, Iter = 72, LP = -29098.6248307068Progress est. = 18.18%, LPchange = 15, Iter = 73, LP = -29093.121714257Progress est. = 18.18%, LPchange = 15, Iter = 74, LP = -29093.121714257Progress est. = 19.36%, LPchange = 13, Iter = 75, LP = -29074.654974731Progress est. = 19.36%, LPchange = 13, Iter = 76, LP = -29074.654974731Progress est. = 20.43%, LPchange = 11, Iter = 77, LP = -29063.4414285943Progress est. = 20.43%, LPchange = 11, Iter = 78, LP = -29063.4414285943Progress est. = 21.12%, LPchange = 10, Iter = 79, LP = -29056.4001743322Progress est. = 21.12%, LPchange = 10, Iter = 80, LP = -29056.4001743322Progress est. = 22.12%, LPchange = 9.2, Iter = 81, LP = -29048.2397922963Progress est. = 22.12%, LPchange = 9.2, Iter = 82, LP = -29048.2397922963Progress est. = 23.22%, LPchange = 8, Iter = 83, LP = -29041.3083204685Progress est. = 23.22%, LPchange = 8, Iter = 84, LP = -29041.3083204685Progress est. = 23.9%, LPchange = 7.4, Iter = 85, LP = -29038.2472940561Progress est. = 23.9%, LPchange = 7.4, Iter = 86, LP = -29038.1236520476Progress est. = 25.23%, LPchange = 6.2, Iter = 87, LP = -29034.0912973555Progress est. = 25.18%, LPchange = 6.3, Iter = 88, LP = -29032.9399038569Progress est. = 25.94%, LPchange = 5.7, Iter = 89, LP = -29030.8134643813Progress est. = 25.83%, LPchange = 5.8, Iter = 90, LP = -29028.4436416755Progress est. = 26.64%, LPchange = 5.2, Iter = 91, LP = -29028.174982535Progress est. = 26.57%, LPchange = 5.3, Iter = 92, LP = -29026.8315486964Progress est. = 27.59%, LPchange = 4.6, Iter = 93, LP = -29025.3351856527Progress est. = 27.72%, LPchange = 4.6, Iter = 94, LP = -29024.8811437858Progress est. = 28.05%, LPchange = 4.4, Iter = 95, LP = -29024.8811437858Progress est. = 28.75%, LPchange = 4, Iter = 96, LP = -29023.569978089Progress est. = 29.43%, LPchange = 3.7, Iter = 97, LP = -29023.569978089Progress est. = 29.74%, LPchange = 3.5, Iter = 98, LP = -29023.569978089Progress est. = 30.55%, LPchange = 3.2, Iter = 99, LP = -29022.7482353549Progress est. = 30.66%, LPchange = 3.1, Iter = 100, LP = -29022.5451488353Progress est. = 31.38%, LPchange = 2.9, Iter = 101, LP = -29022.3619101638Progress est. = 32.3%, LPchange = 2.6, Iter = 102, LP = -29021.9599780928Progress est. = 32.89%, LPchange = 2.4, Iter = 103, LP = -29021.9369051458Progress est. = 32.85%, LPchange = 2.4, Iter = 104, LP = -29021.5441584258Progress est. = 35.18%, LPchange = 1.8, Iter = 105, LP = -29021.2800396299Progress est. = 35.15%, LPchange = 1.8, Iter = 106, LP = -29021.1068901633Progress est. = 36.99%, LPchange = 1.4, Iter = 107, LP = -29020.9971154541Progress est. = 36.97%, LPchange = 1.4, Iter = 108, LP = -29020.8759802369Progress est. = 38.37%, LPchange = 1.2, Iter = 109, LP = -29020.731267989Progress est. = 38.34%, LPchange = 1.2, Iter = 110, LP = -29020.600428225Progress est. = 40.36%, LPchange = 0.93, Iter = 111, LP = -29020.4867027517Progress est. = 40.34%, LPchange = 0.93, Iter = 112, LP = -29020.4155396148Progress est. = 42.58%, LPchange = 0.7, Iter = 113, LP = -29020.3410489487Progress est. = 42.55%, LPchange = 0.7, Iter = 114, LP = -29020.2552523506Progress est. = 43.78%, LPchange = 0.6, Iter = 115, LP = -29020.2188582733Progress est. = 43.83%, LPchange = 0.6, Iter = 116, LP = -29020.2010420232Progress est. = 45.85%, LPchange = 0.46, Iter = 117, LP = -29020.2010420232Progress est. = 46.53%, LPchange = 0.42, Iter = 118, LP = -29020.2010420232Progress est. = 47.98%, LPchange = 0.35, Iter = 119, LP = -29020.2010420232Progress est. = 49.98%, LPchange = 0.27, Iter = 120, LP = -29020.2010420232Progress est. = 50.15%, LPchange = 0.27, Iter = 121, LP = -29020.1072986246Progress est. = 51.6%, LPchange = 0.22, Iter = 122, LP = -29020.1072986246Progress est. = 53.59%, LPchange = 0.17, Iter = 123, LP = -29020.1072986246Progress est. = 54%, LPchange = 0.17, Iter = 124, LP = -29019.913604736Progress est. = 54%, LPchange = 0.17, Iter = 125, LP = -29019.913604736Progress est. = 56.16%, LPchange = 0.13, Iter = 126, LP = -29019.7860699052Progress est. = 56.16%, LPchange = 0.13, Iter = 127, LP = -29019.7860699052Progress est. = 55.94%, LPchange = 0.13, Iter = 128, LP = -29019.6833001257Progress est. = 57.83%, LPchange = 0.1, Iter = 129, LP = -29019.6833001257Progress est. = 58.09%, LPchange = 0.099, Iter = 130, LP = -29019.5814828205Progress est. = 58.6%, LPchange = 0.093, Iter = 131, LP = -29019.5814828205Progress est. = 59.65%, LPchange = 0.081, Iter = 132, LP = -29019.5245687246Progress est. = 59.61%, LPchange = 0.082, Iter = 133, LP = -29019.4905823947Progress est. = 60.92%, LPchange = 0.069, Iter = 134, LP = -29019.4705997676Progress est. = 61.84%, LPchange = 0.062, Iter = 135, LP = -29019.4329138561Progress est. = 62.4%, LPchange = 0.057, Iter = 136, LP = -29019.3843904539Progress est. = 62.75%, LPchange = 0.055, Iter = 137, LP = -29019.3499000432Progress est. = 63.21%, LPchange = 0.052, Iter = 138, LP = -29019.322209426Progress est. = 63.82%, LPchange = 0.048, Iter = 139, LP = -29019.2929128789Progress est. = 64.43%, LPchange = 0.044, Iter = 140, LP = -29019.2684225285Progress est. = 64.99%, LPchange = 0.041, Iter = 141, LP = -29019.245059707Progress est. = 65.32%, LPchange = 0.04, Iter = 142, LP = -29019.2251042003Progress est. = 65.74%, LPchange = 0.038, Iter = 143, LP = -29019.210837731Progress est. = 66.29%, LPchange = 0.035, Iter = 144, LP = -29019.2016981171Progress est. = 66.52%, LPchange = 0.034, Iter = 145, LP = -29019.1955889129Progress est. = 66.65%, LPchange = 0.034, Iter = 146, LP = -29019.1945288245Progress est. = 66.65%, LPchange = 0.034, Iter = 147, LP = -29019.1945288245Progress est. = 66.65%, LPchange = 0.034, Iter = 148, LP = -29019.1945288245Progress est. = 66.54%, LPchange = 0.034, Iter = 149, LP = -29019.1803498298Progress est. = 66.47%, LPchange = 0.034, Iter = 150, LP = -29019.1714165663Progress est. = 67.22%, LPchange = 0.031, Iter = 151, LP = -29019.1698142355Progress est. = 67.22%, LPchange = 0.031, Iter = 152, LP = -29019.1698142355Progress est. = 67.22%, LPchange = 0.031, Iter = 153, LP = -29019.1698142355Progress est. = 69.05%, LPchange = 0.025, Iter = 154, LP = -29019.1698142355Progress est. = 68.86%, LPchange = 0.025, Iter = 155, LP = -29019.1518537764Progress est. = 70.25%, LPchange = 0.021, Iter = 156, LP = -29019.1468486089Progress est. = 70.25%, LPchange = 0.021, Iter = 157, LP = -29019.1468486089Progress est. = 71.64%, LPchange = 0.018, Iter = 158, LP = -29019.1468486089Progress est. = 71.58%, LPchange = 0.018, Iter = 159, LP = -29019.1427801695Progress est. = 73.06%, LPchange = 0.015, Iter = 160, LP = -29019.1329084462Progress est. = 73.05%, LPchange = 0.015, Iter = 161, LP = -29019.1323731569Progress est. = 73.99%, LPchange = 0.013, Iter = 162, LP = -29019.1255342734Progress est. = 74.62%, LPchange = 0.012, Iter = 163, LP = -29019.1222000896Progress est. = 74.95%, LPchange = 0.012, Iter = 164, LP = -29019.1172983951Progress est. = 75.76%, LPchange = 0.011, Iter = 165, LP = -29019.113825688Progress est. = 76.93%, LPchange = 0.0092, Iter = 166, LP = -29019.1089002814Progress est. = 77.86%, LPchange = 0.0082, Iter = 167, LP = -29019.1050456123Progress est. = 78.67%, LPchange = 0.0074, Iter = 168, LP = -29019.1010193292Progress est. = 79.6%, LPchange = 0.0066, Iter = 169, LP = -29019.0963692312Progress est. = 80.52%, LPchange = 0.0058, Iter = 170, LP = -29019.0932784205Progress est. = 81.65%, LPchange = 0.0051, Iter = 171, LP = -29019.0932784205Progress est. = 82.77%, LPchange = 0.0044, Iter = 172, LP = -29019.0932784205Progress est. = 83.68%, LPchange = 0.0039, Iter = 173, LP = -29019.0932784205Progress est. = 83.05%, LPchange = 0.0042, Iter = 174, LP = -29019.0743809382Progress est. = 83.44%, LPchange = 0.004, Iter = 175, LP = -29019.0743809382Progress est. = 83.51%, LPchange = 0.004, Iter = 176, LP = -29019.0743809382Progress est. = 83.12%, LPchange = 0.0042, Iter = 177, LP = -29019.0684084401Progress est. = 83.12%, LPchange = 0.0042, Iter = 178, LP = -29019.0684084401Progress est. = 83.33%, LPchange = 0.0041, Iter = 179, LP = -29019.0574594347Progress est. = 83.92%, LPchange = 0.0038, Iter = 180, LP = -29019.0574594347Progress est. = 83.44%, LPchange = 0.004, Iter = 181, LP = -29019.0487101441Progress est. = 83.3%, LPchange = 0.0041, Iter = 182, LP = -29019.0465112641Progress est. = 83.05%, LPchange = 0.0042, Iter = 183, LP = -29019.0425754932Progress est. = 82.7%, LPchange = 0.0044, Iter = 184, LP = -29019.0368126457Progress est. = 83.74%, LPchange = 0.0039, Iter = 185, LP = -29019.0351594868Progress est. = 84.04%, LPchange = 0.0037, Iter = 186, LP = -29019.0345330459Progress est. = 83.75%, LPchange = 0.0039, Iter = 187, LP = -29019.0303179449Progress est. = 83.65%, LPchange = 0.0039, Iter = 188, LP = -29019.0289460925Progress est. = 83.9%, LPchange = 0.0038, Iter = 189, LP = -29019.028424272Progress est. = 84.57%, LPchange = 0.0035, Iter = 190, LP = -29019.0278798078Progress est. = 84.53%, LPchange = 0.0035, Iter = 191, LP = -29019.0268202649Progress est. = 85%, LPchange = 0.0033, Iter = 192, LP = -29019.0260506062Progress est. = 85.27%, LPchange = 0.0032, Iter = 193, LP = -29019.0260506062Progress est. = 85.69%, LPchange = 0.003, Iter = 194, LP = -29019.0260506062Progress est. = 85.99%, LPchange = 0.0029, Iter = 195, LP = -29019.0260506062Progress est. = 86.45%, LPchange = 0.0028, Iter = 196, LP = -29019.0260506062Progress est. = 86.83%, LPchange = 0.0026, Iter = 197, LP = -29019.0260506062Progress est. = 87.24%, LPchange = 0.0025, Iter = 198, LP = -29019.0260506062Progress est. = 87.75%, LPchange = 0.0023, Iter = 199, LP = -29019.0260506062Progress est. = 88.11%, LPchange = 0.0022, Iter = 200, LP = -29019.0260506062Progress est. = 88.11%, LPchange = 0.0022, Iter = 201, LP = -29019.0260506062Progress est. = 88.11%, LPchange = 0.0022, Iter = 202, LP = -29019.0260506062Progress est. = 88.1%, LPchange = 0.0022, Iter = 203, LP = -29019.0259884974Progress est. = 90.71%, LPchange = 0.0016, Iter = 204, LP = -29019.0259884974Progress est. = 90.69%, LPchange = 0.0016, Iter = 205, LP = -29019.0258648474Progress est. = 90.65%, LPchange = 0.0016, Iter = 206, LP = -29019.0255917766Progress est. = 91.68%, LPchange = 0.0014, Iter = 207, LP = -29019.0255917766Progress est. = 91.68%, LPchange = 0.0014, Iter = 208, LP = -29019.0255917766Progress est. = 93.94%, LPchange = 0.0011, Iter = 209, LP = -29019.025231494Progress est. = 93.94%, LPchange = 0.0011, Iter = 210, LP = -29019.025231494Progress est. = 96.44%, LPchange = 0.00078, Iter = 211, LP = -29019.025201443Progress est. = 97.19%, LPchange = 0.00071, Iter = 212, LP = -29019.0251239136Progress est. = 98.8%, LPchange = 0.00058, Iter = 213, LP = -29019.0251239136Progress est. = 100%, LPchange = 0.00039, Iter = 214, LP = -29019.0250367109 + ## Converged -- lp change within tol(5e-04) ``` r -#some summary stuff: -plot(dat$Ability,(fit$pars$Ability-dat$Ability)^2) #Ability error given Ability +#get normalized parameter estimates +bigirtpars <- normaliseIRT(B = bigirtfit$itemPars$B, A = bigirtfit$itemPars$A, Ability = bigirtfit$pars$Ability) + +#combine parameter outputs +itempars <- rbind( + data.table(Estimator = 'mirt', Parameter = 'A', Estimate=c(mirtpars$A), TrueValue=c(trueitempars$A)), + data.table(Estimator = 'mirt', Parameter = 'B', Estimate=c(mirtpars$B), TrueValue=c(trueitempars$B)), + data.table(Estimator = 'mirt', Parameter = 'C', Estimate=c(mirtItempars$g), TrueValue=c(dat$C)), + data.table(Estimator = 'bigIRT', Parameter = 'A', Estimate=c(bigirtpars$A), TrueValue=c(trueitempars$A)), + data.table(Estimator = 'bigIRT', Parameter = 'B', Estimate=c(bigirtpars$B), TrueValue=c(trueitempars$B)), + data.table(Estimator = 'bigIRT', Parameter = 'C', Estimate=c(bigirtfit$itemPars$C), TrueValue=c(dat$C)) +) + +personpars <- rbind( + data.table(Estimator = 'mirt', Parameter = 'Ability', Estimate=c(mirtpars$Ability), TrueValue=c(dat$Ability)), + data.table(Estimator = 'bigIRT', Parameter = 'Ability', Estimate=c(bigirtpars$Ability), TrueValue=c(dat$Ability)) +) + +pars <- rbind(itempars,personpars) + +#plot parameter estimates +ggplot(pars,aes(x=TrueValue,y=Estimate-TrueValue,color=Estimator)) + + geom_point() + + facet_wrap(~Parameter,scales='free')+ + theme_bw()+ + geom_hline(yintercept=0,linetype='dashed') ``` ![](readme_files/figure-gfm/unnamed-chunk-1-1.png) ``` r -sqrt(mean((fit$pars$Ability-dat$Ability)^2)) #rms error stat +#plot average absolute error +ggplot(pars,aes(x=TrueValue,y=abs(Estimate-TrueValue),colour=Estimator, fill=Estimator))+ + geom_smooth()+ + facet_wrap(~Parameter,scales='free')+ + theme_bw()+ + coord_cartesian(ylim=c(0,NA))+ + geom_hline(yintercept=0,linetype='dashed') +``` -#correlations of estimated vs true -cor(data.frame(True=dat$Ability,Est=fit$pars$Ability)) -cor(data.frame(True=dat$A,Est=fit$pars$A)) -cor(data.frame(True=dat$B,Est=fit$pars$B)) + ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")' + +![](readme_files/figure-gfm/unnamed-chunk-1-2.png) + +``` r +#plot average error +ggplot(pars,aes(x=TrueValue,y=Estimate-TrueValue,colour=Estimator, fill=Estimator))+ + geom_smooth()+ + facet_wrap(~Parameter,scales='free')+ + theme_bw()+geom_hline(yintercept=0,linetype='dashed') ``` + + ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")' + +![](readme_files/figure-gfm/unnamed-chunk-1-3.png) diff --git a/readme_files/figure-gfm/unnamed-chunk-1-1.png b/readme_files/figure-gfm/unnamed-chunk-1-1.png index 78c56d17330df81a870e010172c9e108867d3546..574ddfbafc7e9a4a3dcd6a0eb703645d680cae0a 100644 GIT binary patch literal 11368 zcma)i2UHVXwC*GULPro(K$;@GC{=>MLW5U7H5 z6MsRP0-+0tNR!Z%-i0^*_r81Idhe~b-mG&bYcl(M=gdBB?{A-Qg0Vgas{kti02~Gv z&YJ=N8U_H=4`vMVkCCY&6Y^k9xJ1-N9svLc0FD3-NwI(T>pZZx=eI|$j6ui$tjy8zUo1yltmEIS(yPgVFy&M;k>){}n2XZrP8;0kKGHkmkfX##YA04p+tw?`PQaXV~jx3}$2u z4i>r@75>GFE{rxRY$+_9oGhH&uXF?|9dVVzv6aL7mHYd34|M7V2V0_zTK4w`GwcWJ z9t;j5MUnpohhuSw$ngH~@c#Z}q1)vC{@mQ$rk5)MvELHCKODP1T>1CqWa0i_Bn05E z;Qs!8wlZviXqV;g1&jLtz}fIOqnh2Rw*WvCFgTAV`lK$Tg}*qe@7r`WI;-hATbvYl zM1vt@>Zmvm>IH)&J_Jg=E|_IxZ&g{MR`Tq-`85j!dF56?Om$&($Lradjlk@Y=UcNW zzecYFdafH_@qBaHWRz%&1i)%8|G(v;lc3#T%3T8@+xk`s4EwyeerIwmZ*-II&)aJ#CKN$8BjC4{?s5@3ZA^(B(*DJtqm#@F zB(4p{?9bGQB(}U$0Gp9Q8@54Ps+N?!X1m(mZad-cLHDcJ7(fT1*V;vDhprew(aB>A ziouO)-7aWr zDmiJroXMPlPUdHhJO;DXfB@0-w;iAJge;m|&5oeeNdw!P4Z5070b+O9!}hbT0lJVe zQd{wvqySoTMo?hujnjb4Y+0~)k$+|;%@4>)9-~F>cbgo^%Bb7wE@m1wMYNPJ47zTw zb|!t{fUabl1KVTChg@v{lCC8q%*F;?5eEub|GgKufo~@ z4tA1kj@<+-pK?P5SN@NmuEjvb-@m*zT<`#?9H2&C-VdzboDLhA;u5wD0cfN?E>e$* zM%I*Qb*>~l_>o|)d*&@v@*l@$$smFl!8%1b>(&dkS zOZ+c<2Bl1_{7+nLUXSNvh6?&+0HJ2Q!3jT~HJOE$X(poH7sbJkkjHZjqC1WFaz}9M z8;VpO%}p!SmI90laHyyie2-ujka+v_NGTIckzp)Xy7`2w;1&11&jJu98dhj}W_t#R z;sHEQJ;DsB04jz+?cx|>nT6FpY>RvI80q|`9|A=6qhm>`&NMt9+WB1$&^y8k zx;9-9PARdtl!3N2A*j?6HH29hPXrm}di5FV``xQQ5RY7?U&v$=5c4&wZ;K*396(;0Ixud8?Qz&yayg-u4 zYj(oSn|xoePC^2p+OVMdT*tb2EKk4<`G@Rd!#mGCST`Gcc+A<#2PhbL3Rt#a>#pBs zynsIq=1p1L4{@nii(&`wxz!O_-slMLRgx!{~ zb4IP-O;A=#S<>zeXQZ47H2_^o64XpDE0I2Ix%yp3ec508om~k^3J{uNrgB9qC=wRvdq(cHlbif*b^O32e%&$^v4`%n}v zmkA7vQaIevK<^RhV38UCSFSCe_9&{!&Q&JDe{KWh1OCR{pA;sxU?)`KZ)ZL5xg)S( z8PzyT%(?pe0h%DOo+@r?{^Yj2H$Xh@4(=q4Fp&yJgM={jt}S@}iQmiT>n^7EM9PK1 z4QssR?Dci}Qexf(#eL(WH|+0k(LZtVWGtNOg*DRS6Rx?=_QZAqV^dIxvJE9q z)yeY?pr^UB;5gb~9*dJJq{;o`*a-4g5+Y{cLy7oB4q|nnYz-Nlz(JIMJygXJE-z z?xNqiRl7}*C(1ci9p=qj)h7LcyaKQRvva4r2D%1c2-!*!@8|j=J!$ z>d@!Ute_dXZM-d?7nP7GE)6db!o9}zk)2BhbPyr4<~NXQptimfMuT^+x4M)7)a?F2 zi|hbPg(_xMRbhD;Cz)gyJ-$E94jzxX0-ksa@jD&eqDFf_60Q>`TxPTb<_024(S(^8 z9ZrIb#(6xrVHhmYr9wYXy6Mx%0T|gQT=>S8aq+P!Xw`mMuTd&}Y0L_K2vyeV_Y}yU z(=<{%bA+b1TWR)P5K!rve+cgxYF-ogwX1G+RIT9{gSUCN(fwB@7c2zq*qeACMknSm zOSSuBprI#KcL-aPMd@0BM+MA=S8EPS9UvH4tU{dTGh`;)5CO7|R{{}yQT3Z0Op7YY zuc_k!Jwv-NLyM9pb4g=)qv*^DeL$tQoVX_PCoh;XRH(6U=4lyX<_s`=9s z@x-(yzribV7onSs!YqdZ5~jUmwHnrKP6Ep869Hj~bl|5_x}4xz>GnDKKu-!_VlhM%8FDBXwBrh4W-= zY5+k#M;G9W#<0G^?}Vt66-RCvTHW;EdRwNW(p^IPshDb8u{yq{x99w;DkkBM^*bh2 z#(Rw9O|^w`3GLytl7Bv}D3oup^#`w4^Fg7se8vp6WC zNRKaLrM&ry)cte!-CJ`{!y^n&NF9Y)OL&XpP3pHb0M}5nH($f6!{wa1e~xJRN!olF`G09NHxE1BQu@&^oYCrnYv-J=h2;%_QBkRCKd7*?$3~DAgs@ zmtopmBh&RL>*T~gq}z8-i1p@QJ9LK+Ds8QuH`yZfsP)QU6t3rE$r24|w+f*V|IQHUgx9PD$}CpJHGK;;)mxe~lJYy} zE~!A+;%n|0ccM0-?6d-?dcAS+`j|_EH1>JFN%iz#|3KsYl*gw4T>GU(ml+v?sn{LK zED2W(VTA7qKuk{MHRJrK>Eh-aAuWy};yUs<0&Df`&a~fAg*MNQ-`%`cR!PK?E@+^r z+Kv^_tWtOpzCykMO1SRoWY#HvvhVsGpf~4Q1RpeeQ&m_4xMU8_>O|<+!ATk8om^A`za^XFGqT`~cwZu9eurRHA{T1Q0t=i-<1&rfJFMC$ zg+4?x-~k5S_W6sc)W7;}dQw#(Ii&6Q+@yL#=RFf13|ZT$&a)8mNWN4XB#R+l3V2Oh zf&{FMhvGBR&kP&)?3zCPDu~z_%gIJuBrLdgMGLOB@!4_eO2`F(UNSP|vEW~m7q3V3wCQXS8I%HGj+4^?pa7}y--c&!v5 zZ>J9CfxR3HPq@OoT8}^}7d)z|hnQqBDxqA?XyO}LNT50q=*`aWd;WdSpnLsmVByuy zXMG#v4lOnShhO<%E-F!LU+<#n8F)``Vk2#kUcGrT{+};V6t9DsO%F7HS;qwYD{flk zNMRvB#I1W5&k05;?c|!a0lfg}Di<8%R#zV_Kt53C=D#&$%r`rUeYV926R&_%g;(E& zOE#PPauy{ep5{H{X);4A^E~3)`xBM;f}ah36m+*5T{MDCD`|sXGDRuLTC!3LD-_lSOnf0pk7Z96diQ=U2#FvNPlx9R*_H0d!PR4u3# zNC?v51WtHlRDYF$oh4cMka+AN_aN#FduG zx|>_Hs&~kQPmezDhjB8#)^GoFl41!|Y{UJw+JJ+=BKatfG@3rx+L~5&SSRoIr}&wk z)XVM-awcCyLBYM2Jb4sOnRJzx)^m7&dOdA=qqv8w@wPv^DqZ&J!On9Vr_cX6R`NkE z&~@BJ1u&mD1gs7&XodxkY_W-4%2+hfudemKb4<003x~ripR#d$Z1cn#_2li%mOD39 zNf*^SG#lkPrJlPUxXzJs?NjYF)~4!4L!5()+SD$mBYp-GFI<74q5}kqjFsH@EM@!Er%?rX{pogT zT6VnctHbW5fe?JdmcLx`xy|vzyS%n9_(!Vg#sL2no}=+^MKNYZ8!?7JzggP@TKwj#!KHk(z5UH zZNb*lXU*)mM;b{d;_#VL7rkx$c z6jAY-OY>`TT2ec0LV$A#nr7^hdLO`*t#d<`uTB7gep#gw!q6j&`zJ;`xC#^j`pS-a z#^~<>xlX5&8xkG<4Dhf1O9pJogB5O4qTWV_gU&z~1qFvq~R2h>mX%{Vc@ zb%vOsv~_a2&9Gf&}>EY~b@E-1Ts9KUs>eEySQu;BaOQzIiC)kbq)fvKV+ zjsoXjB=WX8SD!bodA>H9Jl|El`_H+H2LJ(iwfI$U;gTz?)Ym=2C_3Yy@#zIJBkbE~ z;@ROHErFsLG#C%`X-yq}kmNqz>ax3gfK;5g4f}x&fwUWTl_!{;X(Z{eun`po-tvCi z+J3C(n&PUr9wv3S z)1cZX#;Dp~vtKxm1f zNKE-?h;P_#tuL@Nb2YfDzPgKxZ>@5#z5ToFX_X<9`t|!G9w(6MznW?P2h$M~XxuHp z2%66wg##BhSDS>EUjTj+w6_|^%fQ=;0W|k7XLcpj#!ch)P28;qPLP5hqy5R~=l#du z4Y5uyGY%mAT+sL3pxT$p9X1<_9=e+9z0DQ0zja--^8G*ByL8urFS4_rrp3ceF8Mj! zDg9pNY8y8wu~+V1B@G9L8V09sG)t{RahJ2CAujlQD%ZA|Fkz?2gN4j+DvK}PlUSnc zC`XRl)?LJb86)rYbXEakAkqY5o|#W&652myj~aa;32DKB^;TFC*VNWEN*v3V71XUL zN8y(r(&NcKYqQOUKcM2fDyOnMzd`Ho9qg>WUq+qaqetDTi}H#~ec=-puPb_?I6s~& zx)^Y#XlHe@dYmG>+jFO096Ok5o!4;b7gX(@PJzF$`aec9$Y;L^lLBdHSv4&fLJr0G zXj3Ob{OZ28NtRkPmweW7_LI`O*r1;4o<47VV*gCfW1YcPlglSKJCh#EiC7Ysmm&qsXDHw3X&J`J=Z5pQBo#j;ojy5CNw2+(*oOB(&8go*#eCtf z@qboQJIrM+Q)DDHeWgnqGY8lCi;oamm08g$A8J?`6Vk2B+OT{k9ZQ$3!zpEH(WWe) zG%w8n4zd{+4wBbwnH|5)N_Rec?^`;|8+)54l)_^^&a53)d)12KA)CP&?DTfmQk3V8 zwKCsHPJ!5IpVy|>(hbB@=O^=?@#x>WNaP{JgteZfoi)5hIZNqqU;1j#e>(;Y!CXQ%=R7|;~ZJJ>$H|TG}IycCY=t zu{df4m2%k@;8_xHD-!gH zaqQILl&PW-4xm7Woi>-&j}0cyc%u(kGS4~m^GUWKrSG8-o?u&xVk0Kpt6KZZi?#O}P{?PYv1UrUuOsH0*1 zR56%%>^4w9?{VLqU7VQHk&g8&f`T@?$Kq;j9@DBtLa+S@S%reuZvuvE!@pxAmj++- zRmw=KtVo@wwDrNGa;7%DEOf~n%MDvbz&u0cwyB0l`4_ExURtqShU7sYQAK6RVWByqU&_h>-Se)Pq2c77Mv$5ZN#4cCZ^)ooq?SyFDP(|0{` zGw=HQ!~FMu(hcIhQN7ZA$5~C8X=QW94VV6y00j)vo$dDcZeV-&y}n|)OQ8;J=8Y)t zzud7n;CJ_Ijal*t{O{_X_2EYDo+Mx%k}&O9B5DQn;sm|qpfnLG8;n-f+T<*n}0jlAA_gfq~{nij@t5nq1H8ZZNf*b5--Fj3tALh`wWa5WC9EwVv2O6=Y0dY=8N z3q@Z-&A!$J(gxr9J%m1D&CM0zb>G|x2C@-8iK$#?=yfh)RrCuns7e7vxC88u7BErNxxiJH8&}pw%W*Mh_=kAg znhf{_K+cVo>ao6cS9OGSnBa&@Y_)!1sEHp~rud=lH4nsj&M1|wtx}zqI8ldejGM5p zxX8D=>mRYVzhE6Met!Kuu!YeY1(&qbOY3|}%BN=TGn+PHPwf1BXM4sf)t@xi9I9FT zbzPm9RlY!SL%6H*Gl87=Ho!0)-+)4JV@D!i9V)74N?Qmgno~-<*E#XN+ebqZEG`EJ zEufs}lHr1l&qP>3z8n;pzS7x#`QdF|{)dWJGqnr`YRrDkpLFxbbVUF*`CLfQQs;m< z6Cu;uoL#wOfg`)5-mmlY7jbX%6J=Lk&_5MsTPQQZrKwFwO+Zx?rzb!YZtge?9 z#Tp@Gn&hSLl71&)f!=xk4LboyY3cJqlhYzDMes=-2Z+TG82!Zh`M4g zxMaS(j>>jUmVB0kp|UOhl;?By1!p3*4N(>6Gs>EiI(&9twqw_d@ne*N0`IGNKPVmN z?RW?_978!vc$?2xbPoWt+<^vQ=o(A7|I!^+G#L2uj!dC7RJz8JhJl@3Erg3ykYJSc zs>|Qz{Y3R`z8HFB$8MgEj8acbH7U3xU_i}B!hn%D>=_sGcT8Uxf!1{n7vMQ62qgkU zMPzOzqS{8zu)z(?H^Jr19JY+efkQ%2b)f>dRpWgWYExp5{)km|VuT8`#wIYmao$CE zfZEhj|4p#Scnl9xi2dtMDIcUNBS|@u2Z#3P&4mrnXNT^+mh0M;c4(JVxWu^reGEGm|uMRx_tdADpb zf*rJ(3}852KiJPk$dcG0TZw;s=r51Pm2c=xUzgYsGvx=jI>#v0q+g3~oFe5Bxdv>a zV=%9OGcAt$Aj;K8yJoB|e_ydFSru10h4TDT(=Up6kCp0=v8Srn$E$X3q|V-BAl_AO zR`0KIl^J<DekAE;RHg380HT>7c*xt(T7c(39 ziRrpJAszST#UkFAJ>dy%sjW^Hc)~qUgYJGJU0dYyZvd=%V395+l<{r!u4Jx#|EsbZ z4sh_0)s3>E3Z@*sk!zIb-|u8wx8i&Djcap$_D)&9^Rnx<<&=5fZcSnB)9Q@1wsv$V z2BhMHISeORrAE~gaVqx3`s}bl8pefAWFKdljJf^{O@%pnAnM zBgA)|wBi0goFkSbii zJVhnWnT;j#*bTEXVVsW${sb!2IfkC}MJHfj<#Ppq3J2UPRWVx2r6GSIoiX9!&7b2~ zjPv4Wvm?h{r~ts4`}ZsWdiE#9a$EoBlWs-}EL=%>(^%Y8@$*=!O&%~byk6^CD<}g8 zJagj)jczuiWM1+3@o<8HM;75_m4;16EZAAk4kLuEF1d6~elxK$gFaHAKi7;AbWc9O&w0qvY>6 zR}NK6idf)!mLflXJy|2}r7a7}8pV17v%e@q3n{m<>Ek1ApZL~LgU!*Dwwn##MqhEz zwshpeO3k@%oq#M=BEi%cJ9a|dX1&Gz<2r$+<^J&!EnTRBULoZ=S|7UYP$czXw8*<) zPqWj2vpImQPad^UM6{Gh+hM{@5O7XZhUfPfFwVh4=wV+K(EAzNH+=i0cak{op^tIirC8y>9_A;wU*29)f+CIc zgBV%@bE9st0Kuu=(9$h|i3#FC&PNOd-Bez~NW-~q*GMJ=qu%_v=ODsDUb&hqZlrC% z(I!3s%h#W)h?o8NOn}HS9&32v0N`BFA^+$VJ0ZHIOEoVr?~2WzoMmeakHX(jtaABdt$^VRqjkwd#YCx+vuMLf=TC0BXwFh?}sa=#u_|39>X7rfmgQ zt73*qF6Tiny#cZ@A0YF!5KAnZw{+~!aak?)5V&5E$P67~B-UzH?7);qr!r^K+zO#O zmsC6D0o_#ji{w}iK%nPrdxg6!kg)Tm#N{ab{O{h5gZ`~JYzDM9FH*ydUz{ddE%rl9o93`+l1c*qK zzho}PZXS(FyK(o#xg0=6!SB~S8mj#y7;}snh>A?-BGe&=91s&i6EJ54CthL*04rjb zO;uUQ{eAx)2UgwN`4ez8$g4sOV1&~I#CM}tlrQcK{zKbZ8v=MCFXvcZsB*P6#t{kS z1Jswdpml^*%C6OiL?V+`@!0hrl<2b@v-VtCL8s0)Irt>9Kux-KRtGyL2^#j#ob(X0 zx(UC6oC^8+8RcAHjG`jW&t-dXOCiulMuP&--qGin$moMmKZ{hAEKrVt6VU~lw`6ZD z?{!6@SVee|f{|PB)BH(xG)T$P9F)bD?}=J5@oLjqLtF%C3|&aPeMAbLC)a$gnRuEY z1hnT^fO%0ByP1?1ObKaj`W=VGnoT|aHC1B`fi%t|lUzdI)ZqDBPuV=Qau4n!>S+fI zxsLF1*?&`l$Dijn6y6FmF23Li7{2iFvg#H>y%**;d3${N2m4u$6#dr_>syC1Rcb)> z$YJ+7F{XcA3#8ZG^z!oY;gzgX+|%>c-|Hk9|MLZImPea9L)J_;0a6T*!eN|l0tByW zA>U15A!cNum}HRJin)iVmLz;bqV*FeXpTa3f+!wMIB~((?GPLKzkBU@<`1q9pHrTS zG1qs9AZpgp1FPi8NNBHe@xKc(14i<|tbH*9VH19hh9fE~(IqKsuqCoi9C2*RdI?#Y zNCL1YHeZk?9@S{s!DIWAd9i{KdEC40BIFRv131&Fk&daN-}3?l_-|nlf57H`tUQo{ zoYP}O+$0v$bV0ET9p{dsNQgOKt)9`Ddh9sj1PCEbMtm(4oX?1i(MJ=96}=|9oNdyO z@)ffeP!7|(Pe;D4zJg>>?feW#>pPjKK25BFj36>PQs#3jV#cM-&v8`$gV}%c?e$5x zsc&LeBRJXH-8U`I>k{4>MRsmz!JiM(3OUX|)D9}SS&XQxN*HAvQ~j&L91K4%QvUw< zhvXq7m;|pfFBW=>c!;CwKqQiw4{u=c;c{ z91b5gOXI3G2U{nuT$8U36-XfYmA!nPy#UAvqqGAi*Ln*OlXQC1`XDU0f>&pk+IKX9 z6!@~0K!dO!vvnLC$|%}%QB{3@@B&efT0SOf z{IPocJfl=xU`mZwK>kBENFvJiAg^i;j{=aQ;%7wYO)z|?gl5GrpvMdH`3glIHb&SwlU9XiyO`3A%&P@dL|4Q@1)OP$-&rTkl3Dl4Ilf`&5zS*7kWJ;>xGlP(V zLqO>NCHKv2{fWm9to4sgC_4vk_hH!*UGAe{2@6=$e~c1+*RNSjj#mO_X~0LjF-E!()JRH#0qIZsK!8 z%=fJ`{!1&rb{8|s7dzcXX9NQ{2?GBi!ax^dqHE;yvm|ud2Giay*6#-}x1bCgKm()_ z0aTi$+wz19eQ0!kcr%kLILIeB3Y!r6z!aHa*_edNZ_|o6f% z6K&;=+3i1rbM>NuKI4&I(`Px_K7!H+-Z&8GwevI6guFmQ9TUn<|InjzO#Y1TpH)3B zsswuiVq%DM(t=cGe%pyWie*~wFxhsa&5Zp~bRQMV>^(_9{|Ddqm-G97=RZ2W2bVd1 UceXnqM_K^`UE}kw&t4DtABqO>2mk;8 literal 7934 zcmd^kWmJ?=_b=cK3?N+=jdX~Vv?CxP10x~QJ%oe`3?(2#gGdY^Jv4%ZFmx#$0z*j* zC0#>F$8hoezhCc%yVkv*?z7fed!O?==RD{6?Y;JS*4`0NEmaD#dt?Lz1QhCO$~puD zL?{9R!UfQ+n-PeP2=R^ZfNJWg+(-ff83FW8% z0vJpt_NKvNV_`Mdu`&d)H!6(arWcE?iN#*mzzAw!uo`S^&2>%9^))tD28)fwVy|xy z?Dh2xd^2<-ZjhT^$Q6U`&1*^Dsu{Ww5Ky-NYlPjdMK%NkEF$X4U_GyNTn2@8UR_+M zgjnWmRP$8hyZ2HVjTZv$2eK;|KX!!q65Mm2vAW+|aXg@n1SVi=|BPL`@ zOB+NMj5LoV!rj0--EIP_)lx3fW0U3YH3AX2^y>7fKTFEf3!KWwpA+yB74lc zgY2ar`w+ZzfTKVlZ2Ekwqd3+x`so0`UU9lH$!%7V{XfkmE==ixaKv;RgQ7V*uUdgV zr4a?_+(Nu~v4ju|rd2_hab^WCO$vKCkpY@JLXTimH)rSw2~-K;=3FGpE^?q0+nx#2lHM?7{hwti(&Jgg!V$6)lETc~YM5 z8*HXs)1Ap10Dy}$(9UwPUt?!^q6)l|WE1Tio~o^vpS0n(!Qgvmhz3D)?+rt`wBSZ< z7FS$J80rUSQX`|}e(z_i$A+=b(B~q9;PMwI5sPnr+EiMEdR;rDplVyPTCTk=3*?DV zg~cdSo)HV!HztKQ9t3&}hf5D-l)hPe0C+%TybNxb7)4^mV)v zitEz?;D)v-*0Z-BUQKCymj7nQ^Twda|BJ|-EKaDtMWEP73&K=h{YTrFYMNOCn`egE zZA9*HCOR>tCnpr9iAvH7ey7ywlN7jh1R zh=$qzC;X(+f-j!?FRArY({sHPZ!kS!>yh{q#41rk&XV>4JkXwreSX6NTSPq70EPC5 z^QoHM;QtU5mYrgCOW@VqjcFs_WBdU=%P16~*K0sH@zsqlt|0vJ6g2b|y8P=55IK8r z9BOs}j&44#pi= zQq|gYA(otP7RhmE1sLO>nv8U>tmkZy&y8@buWf0V^8w@;-UG&$2t4}{W&`uNbPvH8 zU3Lvjww6)xB5oc2g>3(2YkBAE@in1Y}R#bE4*$=nO^odjuBLiz|4*X3Nu)ltpw_}Jo*^LT) zg)K9_gJ20^2IAI|?q+^Zd?*#@K6T-tV60pdpbnQYE|5vCrK163Z@du8_847ZBV?{;hp z+b0_;Ncg=DZXymP>rfPRth0RmZq}O|?kVeP>Nn?esvrD{)y;caRQjq`kU2 zF~4dh#!{FQfjdPR&<5z$J~zf~jl;6|FqDM_hy9*gNTQw4Wa4`A8uRnl&nfv67CN?t zidy7^4q3Je~#JM zog->JG93VpTD+eGFq@Sv&BRqyw%6>KOr0Ni*IdEQmm)9b$8FS9Aj zH8>vear)jdR|zjkq6}74XJl24UT^2aJ~di?Ci?gH3$UF5891#i{Wwblx{O*|GC|}- zt`+XHE2N-4ru~kl2LPqTwg2`59tL6DT&3Kq?H9cRQE3xoO45NStig>i{6u)0V)6{2 zA~4k7wokEA`i@_Vr)cw<1kUF4qnOJ?Q2+AL#CGYGO=XwXod%LXln1S_0{CK0r8h5p zrKI1&Nsn!v<-rTg%Yom>-W;uIn$~?<;q@nQ-TYpgQ80Ncr@uazEXgb%;xh?peiPvM z>!PCB*l1s6yPaYdnlHJ=$od+_e0Wz^P3b6c1{RxMUUh^2D@I# znO@0>2{@%FTS9mlOf|e)4aE9`qR|@3WB{B7>gTV^m@YdAcO$=EKrAd-$8peX88vPX z|DLS3k5zb^-ZFa>grT}vOG@=&ott^eU{2bL;MNRz?a=5zGbAig^~A`7o04t4O)q@B zChZFu=Kh{oXxq?cQObz;#=|q?K9ffd{2i&-0n?SdIMoufEcUaZzq-j8hXl&K<5=IN}1=2fU)U4{lYau@(1Qn*>L;(@qCo#ZQu-u@o;Hns%Ab4 zd{B^Nr(r{rVB2h=0e&_mOj3FaATJJt=;!jcFly>wPCLfkdJs)Vs(tIo@}0_t6b~iLtJ_yYnxP5_ zs2;E(bUEO*z>v|SdRFu_{DlF_*LY~>glaqrS8H4Y_nM6V?qMgrzD7@vR(Z`J0?*9L z$FUb2Sg>@h$g}_veHjie-fnKT^5|KR?^KFN(Q>^^a;@YPcjS5Lp6AR8IscoQrS#Ge zta9x-tRSoco~uU-1`aS~HA$SD4ZUS$TPNlYtMJLyO9(n%;$dYcsm}V5CGSG*o|0hk z zwZ1l=r;70ag8{zgyt+^tY5Sog0!R}beY7{rah)wWKPqyfeF1)-#NIJIXz_2R;U}pu z9c{X`74s3Q2JhE#={Of`kp>&*fBi}$`6sP{wcqSPs~(1`~%oyeoW*r z$P5x3t(7qLnB-x`;8LJV?LMdVfp&hAq#2x2<5X1-`Z(jm#aHUWG!6p=*w)}r;E<@Q z)rCcOn}VJ>)ED_yDoYa%e^Emzk)>^4a};E);G`^*V7v?7L0knL&I+25w+m=Qq9n(& zUgcUui@F%eQ?dQ`4tUpq^DcWnf8_=pZIFFA%+v5miWvL$ho)F-_uR73>HdLD=2!De z&3JDY@rGY!2o4g{pa*OsqSG7Rp_>1$2&|0K{WX-9o_nycl25!-+m7&7wUWX00!$TS z_tzn&{X?!X6l`nX+RZ zkeu*(4ahW)KNnt1uw!Jm$Kgss%;l(CP}6P+f|YYSq&?sxD6eVlXO@~;(k+CBn$kT{ zfucyk?@a~GC8wW#&+zywb@NxN@g-TzEigqr_x)_a!Gya-_(kmUkW38SDzk^RaqM-enH_M zWr?}^*AnOFMFgCCogM$4Km1hLZc}M^yW9ATDlV9OyCoTS=Su84MSTAve_2H({`|mg z;Kwnt7bcnxWKvyIT30%AXhy4sT>7;l_TM$8I1xLVP$YY$^JjvJoqQDVeI`Y|%%47h zd)qUGYWj)%7FBllzSN0-LxEeA0`j^f=d06!y-No-SfE3n3B&bzWc;#jJyOZYe*`jp zti~X|kL|(_&I;dJ{V1T4+wM&PsXtd^*e#ogK*=V2en{h9A69h$zq^+vh$9Pf@OVqy zRN{IRCyd;Ovf@sTd#YtMVqSB-lz?;}eu*0NjAXedvI$9x>{T+2ME^a`n(&-Iv?h(= zEg~hXVRu6`6fcO!>8Z3dv=&>uiw_^1MY_roX$7JRq7|L!%1t>XcOpxY(KN4in`5e9 z(~DmMNwdH&m0cxBnGO$dxO4R7*TG1vt>yG~3P3N$;m%joXUh8DRtREk-ch^FUNXHs z{+1g9(yI@hALDS)zgz>Z8KF=&&mqd*Z$}PYo+xbupuYbO!&T2Zzls5U9IMCZalSH! z_58no45g6C=W$r|<>7)P(HuUr_qpm#TyBsrCHg8r2GdSY&v*s=sW%i1;e^2hFSoM@ z{;t?)-uJeB|F*Rq>7UVTabMafHV8g6Wttf?)zKKc;GwehXhlw!a#K`dK;Hd{v)qj7 zf_m;;aiRz{+BXSr*6Z4sh0aACj_|BmmML~#zRi5C;a-(oIYZBIjcT?27%S zwa%p1FO7N(=4ORnSL}YBKh7R(F7d)P#BKhVi=sNl*&rytcy6cC8}uTz{!B`|5OS@^ z^(LoO5;Ki?#O&$RzRXf6xz?LeVzH(rROGOln7mrF;wfJMd{%NgH1UUV#FCgQZR+tr z6`GIgTaL@F3F=8omSwzk6o4w44vaWT&}s;?(a3g`wf_3lrc)jK@GS@ezRe7pTuv{-oCp|wxOQ+kgihtdYf6l zp_K`8M)%7fs~#)XY1K9Tnv)Z2v&foCjVUaDonZQ$^1=Pn6{(qnSbFjOt@hBd85gDP zHRO!8f53(`v}UPI;wUeeH7Uh!YWw6Dy)vc-YumZ2|DMPBvbS~A^Ml&al5KXM`x8D| zqxOkSzt3$- zVBAN@m2O>bUmNnxTP;{m4cP^0o3m4T{3fl5mlaZ`zMFp`!Bvmo9zc$fPaB*KQw@_p z=`7$!eQt<+Wo*za_yd4yhWp&og5sKOeFf&k8IV%;E--bd&JO|I4&-h%g;ZvC*tf+9q5l<*xQEg3U^4FuJu%OA?NN}UHvBaYQ2i@!&tN>Z z+Kls0vbvrFgsNk%LK*s4cK}~2Gl3$Seq!ZIp4jZioGG1}M%J$V)d5#n`UUyO5x+cr zo5n^J!%C_`>)oUEqn+ao#v(TZ)tB)I_l$O5H*L`jD z1W;ZNz+4d$C>t~gaW>^*hz{dDldsamvWixqi`aST zj%Bl`ZUYuSO;XeiJbO5e|TH$}1$3thHb zxuS0UrflMpGz0Q+ZS0P?$*|Slt`yw3<5ia?lqc=J_M1*`6xDd~naX@%C%QPE%T`mv zDJbAJS?VeFlv`W$aDblHQuCf2QQWarT;mEyu72Ldy$(%9yYI;hr!vEp93&YQ0<0^W z7RGmJZlQ&g(4U*KYu3!B0zV#OaENlM9k2!|tgmc*K`^U!-AZa{D(c+rui}*DS_a(R8EmdL8#(=UK&s_Owj&LV zt;#t?+?0Tt>5-Iop3XCQzE7%RgsZeREojGfl%QnNqI1$@gC}_dS$gLXt5{sP#BX47 z^i{9qw%G@Bld6jbG3v3PngqPTuFZW?HSDhEHil)S&sU9(>O#MBc-fJ?OxQK}!|Ywq zrnR>D9ZhGGbO4_*H*)ZUb5dzbz;yNYYDvwCww}&kO1(ubc9O*Xon#Q=>wCyXk>i~r z`|;&e^Qsz|vXR3{k9&e7(2JV)7Rh&qlsAoEjD|ivWF7$IwWl_%ywnmDS>7CnK&Trxzt)=vhW+92uRhP`llHB zPXOWD4$MZhb&Ec%)y4U3U!G124IK@gCx=I5asP6>Nw_-KH4uC0kA9=iWGUPBbbcQ# z55>zPIxk0IsG6kOM*1esn|*B5#?h;VVm(Q-3>rQ<>@*TVzb0l_Dv~3%P97%@{*sL+ z2S8MZzrAk>p^_7whs^R$x12S=11ELDr`7SSO`O9|F`0*ADkXWC-m@%qrp+$-2b^=o zSxdx;kW2~hKcRk;w?NZzbH9fr{dIDcRjR88y5!LXB2Ms6J7xu%aHHRW<|K&CPmVt` zw^GXTryfcAOGvnH;9riul5)79)STWrTI-WJ*m?@h@L#2uUZE3RTs<{(__-68_{+r? z9yM*IP3_+Iu8RbxCg)a=)_o*+4#cDv%02fYs-CDm?>vk}rYkwWiPFLm`AG`vv_ z`x>b@bBenAD+4}SIXJU!6RP*TJqAYw5?CiG4bu5JkOrp~8z+4RBLI|!zTe;-4;zQ% z7R}-}78nnyk-e`v)JM5eGxcAnE`cpLY72fgUr3R<fj1&Ka-Ig#9fC%+v4bo&a$$qq=cL*{t#t z53_My+!37fLRocHDBt(`cK1hH^D{~uG2I>Dm6!L<-1{TX3w-naHbKdu9%2`AB~_-! zw?STC14Ko)OJK5SW6gy@@Ld3N5o?mmrnKa?gQy7;U3W`EV%0L&Ghp)pm&Hj@$2iosz-CbfF>=lX7vncLT`>S>Qfvn|Z-2IsWMFGV%;$#4xC;FCYsvhnqg5p^vk7e99B4a%N}S>D5zTV5Y;8pfDn9 z05d%N&y(rwgA!doyQzkhFMt*4f5&uM6gkNY`gL0;>SWi6$@>a9`1|D$2k%7otOd~< zN}qbEjHB|Ntw+kcU^1!#Y_lAq@*h|KvG)JVs+!z~%U)5z&09oZfrV&F_*NHNL)|`) zNV^n#?`KRh>ul?eVtb_pDaeSB-(F;HW$XItZ0Kjt$vnP_b*2A4bCr1yo`> zle|CIVSdhDUT%I}Ye0ucA_Pl~P?Vp}XQ0|aLA&#b=f^{482nIHxpw=d!0!C@mhb#_ zBX^7W@$M{Zxu)AuIsvF5+NCeaSPBJ`ihM6M+P5%Y8xW0v|E$jXxN>!I^+xKcU@y2a zVD?c~R(W`yDBC)BXz0bk%33$8Vp!HWtT8HD3A#Of=ATn-xwp7AVIGL`)7ZN|YNNwL zfN-(8v+DH=>}q=7b+o@lo)x{_w*5rw!qnCrBYbFPSgF+&sAPsv$H5H7hTm_~koQ08WRdH6X$ z-u$G(`|fxncFJJC%R<0jL}-4-@ZXU;4ug(I2TU69+($UvUwseTrmghWnZKiDtcPI% zB~`|Z5|~ NR#(wdE_-Hy_+LIN8-oA< diff --git a/readme_files/figure-gfm/unnamed-chunk-1-2.png b/readme_files/figure-gfm/unnamed-chunk-1-2.png new file mode 100644 index 0000000000000000000000000000000000000000..dc5e29ec31b1c031a953dd96b4e7d3d9b5153751 GIT binary patch literal 11696 zcmY*<2|QHo_xQat!&tIY$P(GJRQAboiELqH=|w16#{QNyOI?IgBulc+P0AoDYX~J| zNro&TDTT3SUz6YUzTfZf|NGz1<@3xr=Q-y*=Q;P>`!^@43}lrS`s_z(9))5OGmy5hvPy(9vVgert;l;oCFm6TLppEOdHG_seLl9b-uk})ww zcIG2HYmpd=%EBO`lzKZ|{T?(NOa&Mv;L*$nde-PVFK+G|VP z8%f$5soE#`S-iK;4`s7YzPGpMtAZ;-P0skYp{+jvu)f>>1K;~mTmaw*V5FyG6_P!l zbNj}Z*1-4DgZ-_ypXAGB<=d+J_8OXvR3vA8_BcJ2l*L=Iy(j_rV4f< zH%SDMqB7xEhdE{d0A5lma1W@|{1 zawWnxLqu@<*2&>N0_Pbhl^L?4?j1|@Yt(e3 z7eVqmzW9x{(xR}W9F||s~&A_PD!2-B(or(=oepL*hwTSt(Hpt{Ng%jW;eU zw$i5j8dmE%RB}=KL}4hV~3tBY|l)G?siCay=S0^;sVq%!20kF-=$6K<(WCMR-qGz2oldu{q91)x}~7?h-BmZ zoSmr8K>Pb7Hgr8j;m(xy4Ps-lhSdCCXqg7OyXxwQd1EJ_S?B*)gh!JigL2#oC9Zq& z!IReGSPYHT%K2jR&|CY;qlDQNaOark|O zqrwx`Ohq0BO1_{h#4_&dT|k3&+Xm4!UhKm0{DGe?xbhm|@m~((0#cC=ne8&v==kAB zctLH_W$Q&fh-aRix)gn`IMxhj#nCQ@3z$Rdu$npOLX+(oXl**>-qG;3t;cl{m3dqbig;#iR?1F|jN%tLK({*gG1Q^!i`$*D zLjWl!<=C_6Gi;vv$9ah-eRk#oPbxo<>yST)r7?jC0(~(Tk0a)c$pK4uUL0{-fCq+Z znry=R5MNw0C|-21@r8H{^Ix|kph~;%Hzy2rA~@zSQnp3bhO-rJ3p&MT<(N0yuTf`o zSoLfN5s)DWSlzNuHNRr9ei9JZmS7}3UWH>td0A;gIiYK7Q8S1#`4gCGeOiTIFh6<0 zE!43T9ku-f zo8$ro)hS8MobZbhzt0MBv`xnn`rH0ub`7&=_p^Ql3^$I zM~Whk=_+3uS#M$q10TM>)Sr?bOu6E67(@7r7-{0&l2tN(#Jx~JlLqTl+mUa>nI*<7 z+@PQo<%yaBlwh+z|4KBWIAh=zGtSCeWP*Y88qxA5h{_*irL>|*jwqQ6<3#uBd};Er zNON3;0r7FS9&}2)#2;ZGNftj=6Xu~h?0BAmdpFcE{GMyi+!L8eCbVX7VsfZwAH;5@ z%S8>S3o`th6}}p-i%XSz=9GDq7H5JJ48b*iS1Snxw7E$Xt+ zge^~~DL{Mz$I(fRe2MWP7%Islj%#Nx!ZzIjy|6K>QXB*2h33KP{v)b^Qz3>qxA>Lh zN`AX9eBdG`VCm>6FWz3PpFYA7e7Pp?^4d`+FQ3#`8jKA5Lv|!hQu$@xA-XA^z?I`& z%>xb~GSjm^fr}Ej|NNT+aX7VChp0bteXx;NwH;&lm|8cT*zFW+;V-_Ys>iT`N1H`mWrPN4ls>Y9YOp%w!ON|UEgoBvDgNQ%EpE~vk21oh!pOA!j zDgej8!`EyGTihed9 zp$A%R+32f6nv5}wsB1F$IO6)cWk%0kfRtNbBpjB3*d!RGHr)FcrPYYIzDxaIIgdlf zblK+miAnH4LfuTl19M!TxAGXL5_E?NT?#O|1Wzg>qGSe&S&XDXPe}vpS@2s(fI8H9 zx$9&ivi+b~wM)QtR8&_Km)c_s;jMNl1cCICjHBZYEUv!vytFKt71zDnLRRl2`uhMm z8K~=K6fTvur*}}9EZI>GLCOIK(a18o5m{d+B!pDHfBsAvnfde&+!9kE$2^ZUW(OTA z;eqg6O84a@nEzZ3JK|F^WGtlydEt#q+_VY+DX%t`+R+1IG2rm}`_J^L0H)Rxs9aR1 zgN8cWg!&%Lq-sdb_V_X;gr1LnXz#SAR|wBw69UasQCE+ zC;zWfdaz>mXc4HqSaTfO3Yo%F(uhQO*{s1)JrO>W-N(p4yrQVIZSg7fUeNhyV7N5p zRMm0~FVr6*i6ND%a-G8JI~Yh|iX^!?k=NLgpb$l^C$2zILNGt9a!R&bk;tTOa+Bxt z>bu&cdqLgzyH5E|%GCHBCYu0#?7zeSpWweH*U2_z`F6Lv*hW1W=;)-1Mbp`)Z~7c{ zX4bHTx|Yl^#HvW2y6DD*Zz=Ia-_5G8;CxS?*ff93b1KSF_8jy?*f(m zI$6(>7K5CXzpjb9G0!&p`gZBSpFr+^k44RsGf|wz z33iR_O6%)J0>%QQu&Gt8d9L?USc*Ii1C)rl#pqGGEpjh|wJ8FmuV=Xe@VRQpFmea1 zph!9}+)ci4DaRB7tT!uyZMx}*r>Wwgod!plfe1GaDbaQhmj{WLsWG5k+pT0LmPhP> zkNu(m`DvU3ghD;tOM%3dChWFF=i&uQZZsBUuVz5DajSv@JGm z7jm*HDo%6%OPzd6nrEaWh8hV@gh8_&n73CXpGObuB0QN46w^nr4UhKWN2js1<1&a& zrII1ce-oBr!l!-^S0BneZ2S%*vv|#U80g~&=7c(JdgmV>YJ3q3K|Fy##SwoahDWxW z4r|K7bF$%4^s9?4do_VO<6o`47;sh=s*ul3Rz^}_Jn1drVQoteKwYt~?!3`5%`D8X zpukQ{AYkk`b!$(be|y+^4cFHpiD5cMItQ5;9=OTMK#4YQgAW?y-_mZZ^ka#HWoB4O z773l?fWXUwAL{xe6U?db!2<6i?NX>DHt0?iXYKr790yc?72UQPUbi5A`P`*sOkITs zZsrCs3wTOc>g85&J=|_i*#B|{5b-=?{QxPbot;zj{kt|v7pU9MIz|7*eX+xM`U;}) zRwPgo`|X}jxGsG~@0Qm_iO4|~DSP^rkb2*AzE%2)&#ulU9=yyPrqu^n>J*Fw@A284;mG?j2;i=>*bgZYk!`&Aci!pNLo>?sd=gI={Ly2LNfd^6J2uDKCD z&2M!61*^;J&ldU%CW*NIhTN)33Q-;$%C(wTo3Un-PZR2bI^R_G|XT&bgncJmo<6hx~BA#ksmj8#_p0VczM>WB&&} z@x9ZkClZ)RF*KV-J83oFTyruN2E|E23qObjo&txdQR~Bp5B9~hm~LrsNJyIVQy&+r z^%e3P;8RNvNX&3PRQ)trr}RqQ7P@oVxZv{W2}}_QfK0m6GwLb9U#tUd>_$#FK%K^< zE4NPP6rBK(Cnw7`3~y69DDGdZR~#}9pLBrgQODQ38mfNl0^i!Lx}s*^;+KzojYhf{V+r#KkIyhr@JzcFK?QZ3^%EQUmR-(VMmm2Yl)2uW| zb)R_G9%r4cNFDaAmdlwGq4*kkbe;(Bo9UT6siXx+s?lhWa! zmPT4~FQM4C`&h9AvH#3b{74ZNm%YY)GL8Sx~33HZAm=3mBvf@YJaxj(y=%<&8@Y}-v=Zok=IPq*tU6f12G zd>N-TrkMZbHsVlzUe9D5iPg|yBV>PwlcK$jS!s1+rHGibQ*91yYz7_Z-`fpdnwyQi zef3F~2rBk3d<%in<&CjgIHTG${i6S5u;R=4(xGE33H_Ps`-~ua7A9JY<|SRtkp=!8 zoh-IhRO6~~;K_){n1uF*2|GDrS9P^Nu+}?zn3fe`ckGqj_0u|Y?uUqilN(kal)25hO$NhMfF12&oBqG?u1*ni4XY*(u#TW~E z?)*rtYL!NskEtB`3??8#IhQ@668m!wsg7P()pY6AJ$l>MumHRIdmwH7v=g)|gJ}8* zI}&k9z84F5oX$PLrs`g^3nTWHT(gUO#(^0=>Mh(QM6|2U|GHeoBy|NDd``3Whs9E1r=Zx8&5V+)J5Q3j9Lfvq??7V|aqKpj!-Rkqub%x2TLGk^a$hwaSCqH&<3gnEtlH(l-k z>9z?uDk`?Xc=Bp){U5~s|B&|o&t71Jh>oM}hfsH8cAs_U&b{K;nb&7olJM zHf|O-B(Dv3YyFTt!Hj)l2oCZO2UNG zP{Nz#IGp7Jtr$8Mv&yUV!KWlf>;so`n$~N_v>3Ted5GNcOrU<`LFyS!CTm7?^Z-je zP00N!yP5x-#<^90Nb6Yvp;KH96CW?F0A(;BsWx_%F5gayAe4t4D_r^-h_>kAJM8S; ze)GHSW21iH6KxJ+AN4A@J$^_zC2l9c)*ZpGsWWlJoa}sjOZ&A$RLGxqYzkc6&~tfdt4EYLry{4ghDa^?vzplJSulER-C$n2Z;< z5~`=#=y6b4(L7@ycGdGyNwiH%Rr#%rQ%olUKZ>zqii!@SnuV$Nm27H9WO=nXjcs^{ zT@#FXNO0$50ImBOi8=%!oi?%9`qgJq{P};nP^km`jGwZ_b@R@Vh~x!5|0K{8QS36(b-{O6HK-33!?sNVV$5i6KrV_hl?ap?jrGb8bcP zn^qHidQb|FWaYU&9X9ElM8-#qV8rUjF0tFuxPZt&VrXpJBQKnc2-2@;bmSH+|B~mI z_#B=tS^UV3&LL)!-CLYeNhrq7-1X{5Ql_ys$qu!_cnI|fG6@_XZ%9>SEsn%Td$PW0c-Im$CcUh zi~x%J(e~?Om6NUpkaY`c8b|2Z*Yf)=jlb+x%8nvM$!Gssx24?uCl;ddcdJB`e1>zE zkr{b3&p?r3xV$v{u9SnQKbR{^1b}|wpahhXZP-%{njQ$!F~(mh*z-m5PYt3jrs)WIAp>=vhd+=;yZ)4 zLXUhHJ7!QWRnEKj8c_nJN0VvkSAIQUrTE@^f`esxh4!Upof%aF=!zzojmU!_yXD--}QOSiM|8Rp2RJeTk zFF09|#7@P4C27$Z)J?S1Qx@3oY&yrmoD(kqaj)dkWUNPt2>KTGr&Effr=%YQNvN2T zU>tW|L>#PAZ-O4_-WgUz+)NO|r3PLXMcC(udOHkmJm2w>HwU07&PniVfXb7#&|K{>hW`bFJuNE*??GFpZjRDgoQS z50U>ae>HMm5o0yU#fi+{4aN20XhY@=vblS&Gx+VVQTa-aY6QjC@t&Ho98y$7Tzqt} zXs|ncX_n!rqvBbm_zm2IPQ7w46nqfUXB%$Fq(kg9Krr;r>C%gbz8cz{M4SaH91PX) z5KDLym};dvtL$$u@NUwLmDE!3l4tLX!^%(Q9X`Ob2^If@wfIbO1bIRnYqt^I!OSSh zM#(V$3l8;=H?^r}RVjI8KBbxL^432Y?j>K;yl5Ua4;@nK^JwAXMPALlYI{JBo-_V} z$IE_1X{d=Afh9=E79grQT<^Vi)UHYA^7qt`Y<7L$U8+!hD?QWDe;H}Ldf|so znw_K=Tj&+Pz45Q}F)GRV62B!zwu6n|(f73#U?{qh$05edxSf{5h^jUXkY{+&h2mc~ z@eEC)oh>*(9(KVq(t30yWldNsZ~+JPU#FR~YGHHU0OE=|CWLzrXgcEfOrI%$BG`RU zOEHHE$Q^ssB2w*jy8Jf%2v7M#X`~D*x0Ovmgx=|WQLAHi=7o(z?fiN=f?C8fHIX4% zV2$s+Mod>S&JuXf24QXsz3N+wCgM0czDP_pn-rtRA&1iE-kXa7!SYl-IU5}HmevDR@1fo)=6!H-umbe{| z@<%Ljea1Fj1$QptxN{w5p#nT5TK3-UBMee!>jomx%p%_Qsg^k|028E1{)Fi>dr)H_ z1m?gF_HWFYiHB_mmPB#f4Wr+n(E>s8sO|G#Op{iyT^mzo;Z0b#!X#XUweQ;S2|n9X zVS;beja15-+|YOr%uC>BG&~lh-*Ue(@^>qXa{mX&iO7*b_rRyw-}-&#|7P3jUW|X@ zitSjb8D|P1e@dKYIPNDj7HyYP5GLX93eQ4PYej2#frvO#JAQWqEiHwfR=USR;W1a( z#(-VAp(@TEJ+RwJWN!D2{n+?qV9sdGu?z$HzSaWC2obcJm)D*W83Uh$)bW{uUZCqr zZHSg7;Q8p&4fz~Z#HWh1>~c|#kH}1>hlV`f-y76xw_9=}jMTa-K{Xe2AcIQ#%?~W_ zj550Y{+jOq<);;=K3D2O237uWM|cx`U&ZWsN|=U9v~6L2GPZ)-{WKI;H6kG80x2&} z*QN=;wU$Lm=<9Jv9A7)Jg*!l@Z+Ow3>nL4mD2f!$x#2Uh^j3Iscd66Ktbx=JF%GI# zP^7sKb#m$I`DoB8ru#-&fxGc33kH&C9r2Ahj>x0ci9EgXAWr#7&ND_#rl*_Ug_rfd zHeHUYX@s+Pc^OfHk)%_2ik~r4&L>>zpR)^K*sy98Hf>i@kFDs<)I*M5%uGgdq8Uyp zP16rdqyJv>zSDsfAo6KV-Fb+v01@`jVK;wp?YjSkPw$Nzhw;?|xh0evUJoL5er*BcQy<$V5ytW7Bdl zDnw%tan1!ZbVfB6$4L~-Rndw_?rn;Lbzi!(Ah8EEq#)(+(1h}=l#_4cV9d*@14K)Z zjh-l%D*=rf*?vC_MSd=)C%%La2?y~ar&)qM4swObGN{mwMEy1xboRIgyRje-FpjyM zS*_7UUq4l3zspID3)W*}YRYAVGX+T3WpTzBGz^n(ACZ-I+GHoPuOB4VT1U4>8vi)X z2&65%e-&JPG!sBE4C{BXedExbKDcPXa`Bh1mXXpWRDEEs2B0)Br)sS?!^2zNpmOVX z17QqP*T;bnE=Vq1g1#@&orR)K@3hTMB(2{CjRhzf-4$@{W<}{CR)h_tbDL6WXJ-mY z>&}O37gKk5`RP8z!A)g%!ZQUZ^1mtYa)e>*m0(jP2O2Xjpwez@%1Bm@fp|UTQGe@C zsPel3F4Zq@cF8Jl`Yakfi_iSI#%aQgNu5eT3nhwtNKJ){<>>-mQ{nCG^k;gPIJEN} zWWTAFVEdOd7^p>2L%oHy=F*vts~HNm(w?&WlRI&JefoY(+MI|w)~ZwTkZJI?Y&z`? z?u((bJLZ1Zai@&UYfD;kA(uBT0H1|Ct#56H?uCn!)jv3+-|e(INdW-!(0(s~Ro+r} z%3mk39UpO?_v{~KqInGL=;y7o0O=#pAJRQ3x2usY-ZfB8=@-cL5`14pI%$ zdXG;XegE=)H~Z$Z!H35(xGRD*muxXHgj<` zru#~Q!ox+-71i+Q-oryVd3ft%$l*Gy;865rkXFR^@(YhAt_~kH8ZY#GCPdJIV&u0y z89Xje`CowaHJ(f_`+9o2*fOtpU^{Xc_42WO>J9_XPR!RlrE7e-;(VgK#MQ#0IBmu8n2gh4zLv!MhMj{D&6N9 z{*ysK=-c5?zpq(srG;<@O!(m^6iZZiZZzNjQlwP6)Tv+d09YJJcl@%ff=YX%`_TuY zz|0cNHC#!jp(#7iS^yHiSN!-BNNnIdJLFLg8=s9Z$)2^ZxP40YA@1X)E~Ks-B0hJ5&NA_7j_*R!Rur)qo$Gwdd)u6I{^()Y*Y$#4eJ65Mf*ABPrsjwNEJnjE~e zNil8;kUc>$qx)>jnM;b>-cR*`%&aNPB7QvaPt8hr6&%_s93fWdR-R}1#!^pon4KLX zbAkyZ-Z#I(_xoi|_a)GW_H6OXazK^gU!DHAK2bbU)2&nc__O~|uR2YA(s=LpBLlx>T7|j(PST4Mh#UQ0{5pJEPlX^VF5MG*UoqZbg z3>&Zw)8K0PWJNiOynX{_Zu~;gZ7M-p*uz&WfG#`8M8jUG|B5Lt^^hKPtZrBS1KPH5 ziqhFuutMqIgx!>pV~M7DhAwnJlS10F!LHPj@}rnc{Ho>EqD>cgX6wKzl>_K%jtw;} z@U^?eOdP3Jy^b9-al=y>d!%4HUS*dXxy>#wnNEG4?D^r`PtMic+oF;Zba;n0NL zu+c5@ z9&hgQkg|%S1IRO*<)#($`3**fNgi1(py^l(7p);taRJweC)~E=`@*@nDCUZjx7mL$ zeG{_sWwAb;H9y#!Hy!MCJ^ct015iHDUV>*?krIuevuskrignD=V-+`m@G!k`|BH|{ zVfN9+q!-`k)mr676SRm=Z}silU)`>!pEycDih&$qS$?sNQ3*NSZvh_qxcpzE*GrP} z>l~xtgQ3tGoN-NwSX7AlT*U6VeeAw@%~SS6L>v2(>HIrwHN(fZoYy*8tuH2_p-n}h zusC`--FAi|ip1W@U$!o6XLZD>Ah{Fc5H(P?bR z#nM~-(Or{+j-c*j3G~J13!bqyR3=9Xk9yjWm6srH}h*%C{AxDEU zUrnyfde^y-T{NFZvUi1pm5Y^rGp&3O_>+ zT1L};^WnSeS=aQUa5162C)e*>0=xV*$w&6J%Le+{(Y_(vw`1v5uF}(n;cx&NRY$u6v8ZnOELQ0n7yu3>|%}MWwtLir%x}?lRnc|%GYTll9?a7;eSna6IcTbr8^@)523&# z$mn3d*SrWkdGKGQqybM6J;3m3IRoxy*rX#o`Do-=?cAf1&5!l=xr5LxHz9b2Usv|E zewpSW7i>kXsCSrtdOa87@E=11eTMEFtR(UERTi&YwXpO}({)xj$mAbW=HpWZHM@j$ z=b0%9AM=hOR&t2S+AXKce!{U@=1OqT*?*|S{!4|Yl#(s))(b_+4r6pvSn8t*9!aa1ClN3=lm!w~a+ zX7ayHh)bYrbl<(f*|lS2HF?$8gA%~qRDy+~LHnuQhU#an{z2~b4I|o_o|lY&k>tl5 z#zMZm3RTm`r`kVTIj#|WmibYuQ4Y<+v@&=*!p}A#?t`LcSp8A1CR2@t4>b?(n}_`# zU9$>+hVuce`icdgO7iVrT@+4uKmaLJk;WsM(hDyExS5( zMPKEgH2ongxh4~jHrf@8Z9M1?xkY&YPTwC*6Q4S?^t*cdo-NOk_e)mkC%U^2Jw*D{ z3N8EL|4_zzY6biv2{rcAvv#=t%&FrZ)Qe}Ly<6M1Zco5Huv*d}rS}y}0U^{WzrYF0 z*EtJt>$Xk3@_*=zCdyLF}5(r zr!mAI`DWIq$8)Za+1vb>epe8}HL59fsO#Mp_mY+OklB``z=b*$*_-^*fufu$Qu`Kd zwybIRIR9%sNOU!NV64dG##OQJhMfPe5mXnlKw|!PlZs}>zdAlC4)t1Fvha#nU)wo4 zI)sUVXIvu}XfbTeP18R2yeuOA0hFH9+-Z#y8qC{f0{e>}Ee?h${62Ke=xa(IIi&mD zZl26{vV+1bxYN9Twt=3XC|3FHvGs?EO{$Kb z#{lI{fJTycR$KPOeKuF@&o7t=ahq&>)`ix3fzZ$=9?!^)>xX+|B?EiLtrn3vek%bH zB=icLJ2tE^qIbT7p55b!1%_?!ZtCLCweneRVwW;S3Z0A=e(H_(HtqUjtX>A!8f59j zP3(+)tu-w@_gqJx_&>V$893&7o-!Y^jHR&25TT<;$rRy`p>k5@L#D_Ohfvt&6QMF?N|Or7 zkwagkB61K#=92&R{r-OU-sidZ?q@rXz1F)v>s{|!@4Mgi+3NSj5;^#M9U~*4VSyIJVPv?+UWDO}7oIu}y#D z`7es6ajd7Oh$nItkNW&iXFfjvV)+AVXWo9< zoj)^M>?u<0X!Fre=D&W_q{T2Wa*&XdcgO9^Y-=-R+Du z?i?L`lWzNFcXzbdbF?#dbQEcdl#PyO8XzL$yW``#yVEs6)4RJ13kyFxS~o)j|Geo& zfOg-c?~Z5gjyLa#OxNu0;UO@4;=8-ML3)@+h{;(l9z7ER0Gu!O{=nD4)bjwa53n*b zakzGOvG7K*xaF1CPU&Uc`!H7g`|h*HdK_fFF-^94en2{l<$CBVmOSBL-tp%aV=qk` z8cv*+So=^_vrSzXo|tf)*ol^kxT>v)Y1f`;FLJ*?pE;g79&MQ@$@aSkxOnQ&>xlN#RjS!#ZlEUjr*kGHf~y6JbB&} zV&g?s|9@8K?yMg&es?119C;*TvQwWiJMiRh>QdLBIxL8)o*ViRrG2b50dzJ_6s0BI z99){sGld)tHwTTQH?qaCx3R=u!k>PB%25d>39f$L{nL9P#(w9Q*jK}S7zf9-)k6*Y zA6cUs!`OdBE%QB$uK)Xc^$wS@#KRO5$T7;LldLg=f2<1KdHr->)Dg^T%jccx&lUmY zvT%xIN7PJj&WBe)e#6Q*+WLNipf_Y0gp!6+PXAUuDdElyFM3h`M*M8XFkWO>E}F@} zDMry{+qdqeYwU@4Y*Zv!i>*%jx{HMfbAmZ*U!;hOMSl23aj5+@VQs4Gr`J2i^%nQv z<`R0jh>09HVZEqzU$5LmvS7Itu)Qo~QR|A7>##tUf+*)qZh|Y0XkrRo$>CtUkIyi3 z;y?H-1w$9yycQj6L}v0e@?Z#KWnnB<1`fxgoC*6-)k(P^6<`Qy#6 znC;By$gRKi3(7?t`wzS2auW<2aZqt&4SIDjFyQ)FJegD0oy|GQ6bjq}3l>`)iy%oJ|rmGaP>(V%~%O1Bn~gMFQOO=eHTgDm>V`XReBT2?%F!E$h) zsJ#x#`N9Eca^r;R^^vE-8Qa_T|KQ9aLw*)~ahB^SCAjjjySdsH@$ab%W+iesY7$cu z{4_#oxi6g#`B|i)O<-SSlAPu$uJKw^;LR6#59h4AuDTCo@T0zWvMYZYkK4Y zA5!*{+#CFFnNnZ+Sy|D?e+0l1801oc z8q*vr|AEWlLG;q$g(v=P0MB_JJ{UfMW4QH-&%c#LEW$>p&zb6G5_NmEvmlkHhL?AE z6Iqinx>vtgal zWPv3;S}=JI7}mT51%xMhz1ER&TQHGIkwUOmWMcUnQLQI=ZxY04er@B-PRqW&cTYS` z=Iq?hiRv{@M2pUO8Yw}bMf;b)@jI4~Aq7j|D;Yvl*p}pCZ#EAfMl5!iEQ%|XfG@sWx@}WnLB`I6>O758K1IqEbjba<&v0BgSEW^*9jujirp^8!MHYY+$xJBa zjv_43E2kMAAjc1SuNWv9|FLvnv)n1clGp-{g#H$w4EN zzdTE0E=tF}lG;S2#jW@%idg-ri zcme9IK|h03PG-EHu;B=kIS--oJT!M53nqvZRyHtTGVd$KLxZ3Uvep`kdXwS2Db7TE zcqi#MP;h|WHoiD-IokJh5m&$p3nZ|X$-t9Ea<8)KulH>2_%TPo~<YK*PP%zX|bV#J1m*%WlhNC)uA_vj!PKq(Oh%93I-!P5lI_E3*7F!{?x2cnx!9lcT zTn6QVJ<*>9a_#ZA{`_>b*{Zyo{|IeOua*Rxlt0EjS3A+2F1Rr7fP&3D`H>=CngQ)5 z@!{_&0Th1ps-q#w-;oyJR$b!8#MkLO8~6TTO

*zvO?LJ`*vBHI1a3tPYlwx;hUdPMd*6%#&>t*^`Ik;Rxe1fi+@_m>_k4pt8 z$G`QY^OX8xT>l9Ey0da-cirywVOdUzWrpY>#?=(xmpWjePhNBXJ^ct>o!T!lz|0bN zmX}%ZJDJH40+aT9lcLhKRPAQip}*&In!v=y>5DvMQ$b3#Q*hKCDMNC1vX2?`Iyw7Y;tqGPjxK#7_!DUHxYz@;~sG+~h~l?m;|u@Ud1b zUjz15D=?NfxP*DNQq0jDT0&a1kDU&K}=WSy)RA6gZ^zW<~<6U@ZHhR@bJ zKJU^`ar@@eOdd&opWE-Yq|q(ZL4qzEJX^_87L_}i@|KmBNDK5e+dy2C_60grnG60Y zXOlTcwpco*S9EH1|HG(%LJy@d(v^)4VH|Fyoj6oRG`Y|uYd5|8M2v*G&v>s&{rF4s zYQaW?R2F+)T_9fxj#wi)xajq%gc;|`3TwveLF163l|jXGe+3Mk*QGu4fwp~G;9x6y zz8iLeugmcaDRBGng90ycwX;f(1tu&XXW4Mj4ludK+GRCvOF~FhQy}wQIay^*17m1N zcoj$TRPd_SN+I`7#m6PTW2Y4flbE-UXGt#fEpij=RniC#OD&-?bYUJ#X63;3tWnhV zovsgIH5Xvx6l?cZ~MI-=LVc z^UYL0nV24OK^o)T)KMfqVhZ;P{b#>uTNdT57~ztcgMZ2h*I0hYDEV@6tn1CYML47U z_@YZr8I0d5Z})qKljq|{EE=0I>392>C^4%KubWygWNTv!-qvJD+(%>0hR^)zb|(0} zm1L+syZH`7@an}Qrh=z7Fr<^({(0T*%JBOIV2!ov``C=ej`$`X!R*t}%#6v~QV)l)&Qmkg#MK0DNs*td}9`r*> zVle8I4@fm-czMcDkO{72-Ui1sS(&s<2XZ}41WXe^-=|3CVu%Sc1e?^HjEofeKnEF0 zf+&(D{ARgd)F=d>GV6&00*|}>6I!+Y%mN9sY6K7Z{*+)($~B-@_Y#D3f7OPW#T<%Y z|JZ!fpep~%aPN70@}t+?I9FZw?^q;aOQTviN{xQOK*{+ta5G z1hEp`vIMYqVEw<;hlWJFcrW5zGV9GOENWh+QHM>rI9=|pj(+oN+v-y;pZdPXPp6_t z-U^M6mE8>7(~Hc8h|Yo(Y>IE#s2b7N@;1?=r`bqp6e3!YPsB(J#u>#_eGN$=TgS8B zEHJcsBLjNduKSFfRMgRjtC*SfkGt=5{BF7UJxV(WicDI`fDIizJk049LnbCY_jgISQV>F2A*X+L8Qjxmienqk0 z(4e%yN5c(A7cqy{2sdb=gC{^;;v@2Wa@{scXB^iJc)pMz74q`^wmy5CQkG&|;rORP z?8DuA3hY#CM(+U*YrzX9XZe&z30a@xcrQm2dirO}K}T-aK-&iRd@EU$rJsXZU&(>D zNAf3mHYNw@Mh;WxH7gy!3dz$@o>(w8oRGG5*a|*B^CQfbt;P%nms87A0laNf9!Z>^k^nOB4iT@_a+-YT3r{-OK zl7pSuVJy%ac|j-+d}CT&UAbIokn_|7s+c(C!%l^W-}D$PWxq)If)r{puMq4I*^(3B z+QCC7iB&9~x(qn$n)dXRRvJ9B?*gE$AeE42cp{`H(43A&*VuN|~{9;DL1MEU8P z#X<9p5zF*%CTHzqQ(p4f9}HQ6UDSvY90MHK&7d{fAh;t2&-HMv z$dLZw0(IFy35!oj_{>zPKbXBI@9;}bXRmF=?AjJZcYh%}Bpf_%Y-VyG@!QPruDY5{ zE6M2`JjQ`b7pz`?ly=eJx#;KyRfuyXoFq05>ml?VZ@|%>yA@ERhJ>$ZCM}p(n@viNYZ(Mt#lTpFfN78vt{cx!lyX&*ZM6?0kBg=Ec zcUJYWtNz0+2P+z&s`|*6fty#P8k@yRO`BiL>ruy0?k&o>)~I?0 zavp;6Cp&tcZ0^in0XVz_Lf5V*%w4`lzselOBSNd_`AZSpI(paf4@*reB5AT1eYYWy zUj!!6scGEo#G6yhZz^-N&M28Hx~`M4#|q`XVYtBU@2@JB`T%H$%l>Fp66T@d)iYr{yBq5 zh#L;gQbrBjnP@xCseU1lRvi^SH%pSebG>guGus%x;dq#JxOj@w9xAao@ll!;Pd_T^ zk;jrSe|+UHT*SHd;Bkm{p-|M0zO>)jVxvJ5166Ngeq`{#Pp#j#9=8+ovg!^zIjwxI z*+=}|xW=nv7lWQih4UR3ZI)x6gDEB6E1iZf`o3GfUn(`RZkyWl%$JmZDk0x}My@2{ zrNl9}CsNmTH|k>u<0k(%@)xa$tj2a?=6}e87ox=1^s`^O{;C$mHwyX{wS3tpk*8H^ zvvgG0*(W$fMSm-$jPI)3!qi#LLkqEhwYk7e(xS+^5=8nUvO2 zWK-o+J{V99@4s0|Zmib6msGyvnlvoaBdF8yHZ*uZT+m_MRp;mpC-JlQ=bJdIJ)oRi zdw!Uww&pa{-pLozWHK`$c7zhzJ=k6Cwg9s-Q4Z-3DKlp3JvFn%%dekE%VMs?{$2dq_r}V8)@W5lOh+D`I^G?(qC^=nym&_!t zkZSr&c{Hw+fRL8~d*tN|uo`bNqP}^)c%JX3`xf=ml~an^JOqd06%-{E)l2kl_%80E zPHH)_c5mmeO*W3PG51P_|7G>bzsxaF^=QhD8DxK?56h@2*0*`@)ub_Tge^tXNz(C8 zPe1ER-ap_!@FVi;Obm)RzyGCV z3FgH+b2-+*=kcvBka2C%8zJK!do6@*ZK&TZ-Ck~z^2N6~K*rIv;c`E`W%#2N&YrfX z>>tfARh{8uLu-Z8L9BmY=*7M#2q_A`8`d7%a~RtBM^Ps2+-{h zOvCA>s*FH7dF01=QVt6P6|wx}`6Zf+q*d%!_(+S%T)pWye9%!DKA3|ee(jxJC7MD} z2+2JEI*thlA6x8q_;!G3huFb}Ea^W*esXQYy~gNWScvn+knvxY4LgrZNqE)8zI)m%sO9(i+2u!j`ZCjmuA;-2NoS(WtG;mi8~nzpQ57&0C6YK@iG^0pV{P|x z7_Mx{;3nAO@Iw#C{s?qD+20h(w_u`^)_G9Xet!LzkF-aghHSsWV+YM}RE0*^<~nM= zLDZ__0LCFTtN6}f>2yfTugsr%Z;urTu~DN$u0=O)-cli z{YR*>y%r$)!O9d$z^e{d^3h9NMGV@qlX=%|wh3H?a>%X;f0Wp|8fL~4p%8U6==yem zs9nVeUq}waZ`SCteYiBvWE3dP9l%)OhA}gu+fjIy8cLIbb?{i|Z zBzrE&YnSc;bg$YZo}3=h|EP9RQMvf))B|33Cigkz`(eX?uJ*r;&6s$pYy>&@dzHV5 zfULnuwq7TZqYE7ImZd*x_lt=X^^j9O-np?8AnKxkPz8ARs&3=Kgu3DDWv_CLpl#*Z zt@bCzp~WZj=RXRw6e_`y6)L-CPwjYJ2I_w~yxTF`e>fY=`Do?IB_T(u+}lMz?Y#Ct ztu3*Nkz|2bJJDAQmV*{b6u7DgJApvF1$Vh&e&%$q?As4a0z}OAu?tr`NaD!!d0taM zTTeWmT+PX_LyyA>rZ0dd!~8GPlHHc*UZ=K|<1e*1tF_ury)F#@;?wvqE=3RTsxLH> zD&&?BN0mZ-;^MLkdV7v>iQsR*@fCJ%IB7~!!>o4q6YXsV&L8lbm2*vHc6i#zmN@f^ z^u?p-wHMn*)y^G3iz2J3n{^dG&_SK?_;Um0TuPE_v$vN9Q$ESR9;X&}EG(4)94sxL zp1&dMCHs@(keu0NzM#Bs`sJ+Vne`!7h7SqK(SI_Qxb`q$G71>etVSEiXLLO9eV?=~ z;)9wI{kgzHV5A+g0GTkPV7f^5->5H50zBSKWP?|~w^4Ner)e#jSb}z`N%l)~vz{+O zA8cMIhGt=)#TQYbs5?FG7cPg4{`v8O#D+$Vv09L3eRPbGMNx^uA1?lY5?1ofOY$op|b>bI?Px!sFhxJlZXACOMQWT1x~G|3EUr>l=?hS;Z-b zU+jI!QHjBOb{cp&Zt!T;y5IenYHIC%$AGDNL~86V7>O^tJl`& zwP4oN%?QVy*_MLij2}m4J}<7l`m+wwuSN~o379+le)I?HX3kC*>v)X>m=oq)BriQl zwwF<=<>Jw*N1IDMNyCXgk6=gsT8j_L+pgnhS2{XW?Nb2y7sH?{(eckH3ulE1XtDy7 zaW0IvfQ#Tkau5-osRg?F`gkjF5PA&WewvZQf_Y2<1>8!x&qF6dIo)cKbJ6(y`>R3G z&ZxDY|AXf4HO^^ZzHxcAMed0i&W7WYLBw!kN=vQ?ftYnq)}if6#($6;pT6#cMdhB6 zd+U#3S8cd)nhqWLwD08S@sA=-2K9%9rN<@SK<0#LD{QtvVNL*A6RJQd$NHg8pQi`~ z$>|uwX1a&#j7N3FcGWc(q zjUp9KffsGF}%HoMR2NW!lwsR-3iujYpI zI}XktIc-^0hV{^W>bmS>dJN?tPnxZpMYv}aSc^Gro~C&)T?$@XgOcmUyEyV_rGR75 zxt6OVzYrEVmfX7uqVxST(FC_Pd;*~)p@e6ClK-$CTyVUhA}#wg!Px#ANL>#Y2*~n$ zPF}fi(J`kzEr+%O%s6#@9q5h6qx^HYU_CyfB_ESRC(~xnjk$4^bIHA~TRTr6NVh8D zm=h1!)XRc5%S(mMYF^~1YyKD%j41fnTpOM$*I*d)W$uq$UUHPyv9>lbG|Sh7qX8y> z)$Mvsb@^QTvG*)#Ob#x8R;C$iZXotl!gsZ`W7Z};_6lrN4e}FmLESS)$49GLCA_e7 zy;Gz}+}BCud(!F!lRFnaoJxo{E_AEw7WJ+shoc!->Wrw1qjzB^L=7-0m}IGky}EiFGMs^aHF z2p3m|R(?*grDSLsLc>S5z1A_5pC8pOC$@|}i`RP>9X%>G2EPDVj}M%!oGfKaN!QYa z){geqVo}u!`n?i$yw3HEuViPIrq!v0ZjI}pzk4IbsN+qXXPFHpIGh!0K4zrS%9M2*$beK67=locC84wc1)6olSdiw zm)3c0lf~O@K`o1-|88Ow@WL0%&7OyYbQ?yvmtrFT1wMJJ?o{x%-h4P7&^1#=5oAkH za1gqPjW8t0RcA77QfEF$0(8xl(FnaTRdB+0eH73Q&2{m*_Tq<@fdJI6j6wG5AR(JW z4NIUV`7OWTn_kg7+TnuHQENSQob^QsJzjZQq4 z5H3>Qq3ti(+5~@ior^4_Q}ffZ>}m!J|~oywrMIElmC@3%My_d(si?;638vHz5nl|98CNC1)k zy!op6MA?a@Nlh6ZAA*uz>=5iG%>iXq?)DhV&FCz0Z+xAELL7p@;q#Cv zsjn~i+QqFU``witu7n>|Y(z`O_U2Dv8UFh_)}X`qxCx{K5_a)p99-(3Pa0LX{beMJ z3~Fq&n+DkQq*6G!GCr(`%6oMqJBX+6!(>a30Td1H_G=4Ao~76g#aZY4ab8qZ(rdI~ zvwrEgD_XK-_PEDBL zY1Y>i;1^0>!Fo16irK2M;?mxKC0?(JE{vNH;-H4vwh+y41Jno5_rbxzwYs`ub=Xzj zeW3eJK}qlWrEbt#{CAZG6XCi7JMkFttE#6?NmhGs7};4*f{3JPRV>6KL^kRTn(LQ7 zVF|)`vv72cDL1YB=QBWF2F_tgc)ifaOgE<%4;i9F#f-$^n>0?sb(TT^9;r>};5LW& zo?9Y(5rZFSl5Gp&L%!9=bK&q&^&iXJz{W!yF3Ikyu%79Q!X+MB=1*x*o)?B$bkV)8 z+!{f;e_7bz+h`Ld9i0gFu(e>`3U(B)M&n6TW`k%eaf zRvv#0<@l~qNnZ0%*A{SPNszx(fP;j;KPQ`zOsd9_Q23tl%sNEceplSx zK5-}VCl9T0+77B%voa{LfGP~K|{oXJ}N;HRgn7; zO;9Pxb&Fhx!>{B^5z+|BMhSyksDj!W#_Cz}=+UtE6nYM2t4vv=Jc;{9(eTdHh4b`gTlD0h%{yH}6HQ38ih=%P< zO_Hj#kMT=K?1hz>b${czYer;V=CW{Y#~Nf`9h+r|8e5y5MNuzIR;lOPV@U8cIC@<{ zN0=7icE&ED{`&!KH6EVJ=E%IfMvd_rNU>SFuo;^)`S7bTpq?X-KEF|*@EE8-w(-tM z)+M(TrTb0*f%>PvjSFOsg4G4I{_lXi+rs?Eo2**`C`R_TRrSf2a#|OlaqH+%^7j22 zZ7pY2c#nNKO9*S(&IhZxdu&~@bsfl^m3w4%&T@N177JwtdJFF7y}$Gb@g=7J8E9j^IKes9Z)^Dw#R!;H5- zTXH81W&K4v+YqP~5m@QrB$lPhH<~<)&f+#ugrYtS$J)ocu%ao!PLQrKB+qB9vb7co zuDLnsi$HU`A(G4$Cu|wp2kv+3%0uXp;T$eJRRo2MP=r&P#T1%&RsmM$(ZbN8kLox< z3AwOq>tYpILP$DEoqKy8L)S+V(+$=VQ|&;CvuLd;3Gr2Ux56YJ;%yZNJW>uk;xs{)=huPIO!+fE~pB~}& ze4z~}^S0@M%Hi@ZAKN5#exLQ^Z-j46j)DSDg549tlTnk!3T#B$^W)NLh<6-Tj{pfn zn%N(j)TZ{;(~+KRuH7SFuE+FXFZGhgkVwL8Z6pl;;l{2o7{u*Bnlp}g@tPoXgBgCH zz)77&V2ltLagZZTO9g7+80Szs7nt>PQ3hB$ldo!~E=E`nNf#uwyRcawKBd5RY0u*y zj%tV!W_!;|EAJ_w_j|Dvydq(x@^a7FEq%KJAPJT1glowAZH4i@g>8@Ljh<`;n54-$34A)qP+lHAR+zoI_l~Nc<-(tlhFvLf;7k&} ztYjKMIWDavR7*O_Z>FRJ-k5r*H+K zaXTfBynZ_p`NAa48vm|djg?@(H^ygxfc}Hd&1c8CY3CBi+wn*Y{s;93p|JmY>xbP6 zw3g1^Gv9Sk`FqVD#`^Oyxb)v=Kv}AkfNeJ=J)L9Imv?c|r=xo%$Vi;a8u$8$t;89E&S0 zcbQ5cuNQR4;jMe+-O#8M(+W1wTtE%O@7Tyqa2Pa&?9DH5L#+d?jpyu$o`YAS%AenJ z_8rKaI;r8@KQlElI4EHma8-TE#wk)=BSuMeUBvYWpuBnaGgq|H;fK?0Ms|}9i$2EE z-BoUl1#D}%rH6RRXH3m9b~jzBUvJ872dkes0`%@gO>=~x_)_Z@WCiTv`=d=$j&~Au zZ!7YrJ?_%TejP(nUp;BT0pd6Erp3OQ$qBlAc&*&=u4MF>JiJLyJrzo?WnJqb{j9rRL=?^M1$2gNWt5(=QVkFn9#2ZH9xUX`vIyO0!v4Q zXQRXI)MOJ@M(=E#=Ggr`u#k*`UH=y#@Iuh`*T$tWcdsjv6G_JPf1|80l-DL}(uj&U zTHQCJn6(%`slTIX$2R!cQko*VkZw_cVHZ7gq_}fSdn~Cc(EhUe8G