MATLAB Examples

## Statistics and Probability Letters manuscript

Tables/Examples/Figures/Data
References:
[1] V. Witkovský, G. Wimmer, T. Duby. Logarithmic Lambert WxF
random variables for family of chi-squared distributions and their
applications. Submitted for publication, 2014.

     [2] R.M. Corless, G.H. Gonnet, D.E.G. Hare, G.J. Jeffery, and D.E.
Knuth: On the Lambert W Function. Advances in Computational
Mathematics, vol. 5, 1996.
     [3] P. Getreuer: The Matlab function LAMBERTW.
http://www.math.ucla.edu/~getreuer/lambertw.html
% witkovsky@savba.sk (Viktor Witkovsky) % 28-Aug-2014 09:35:12 

## Description

In order to illustrate some of the numerical calculations required for testing hypothesis on canonical variance components based on the LRT statistic (37), as presented in Example 3, let us consider the following unbalanced one-way random effects ANOVA model, as a special case of a normal linear model with two variance components: Y_ij = mu + b_i + e_ij, i = 1,...,G, j = 1,...,ni, where Y = (Y_11,...,Y_GnG)' is n-dimensional vector of measurements, n = sum_{i=1}^G n_i, mu represents the common mean, b = (b_1,...,b_G)' is a vector of random effects, b ~ N(0,sig^2_1*I_G), and e = (e_11,...,e_GnG)' is n-dimensional vector of measurement errors, e ~ N(0,sig^2_2*I_n).

In particular, for G = 10, and n_1 = 2, n_2 = 4, n_3 = 6,..., n_10 = 20, with n = 110, by spectral decomposition of the matrix W = B'*I_n*B = sum_{i=1}^r rho_i*D_i, we get r = 10 different eigenvalues rho_i with their multiplicities nu_i, i = 1,...,r.

For the true values of the parameters: mu = 0, sig^2_1 = 0.1, and sig^2_2 = 1, we have generated the n-dimensional vector of observations Y with the observed values of the sufficient statitics U = (U_1,...,U_r)'. The true values of the canonical variance components, vartheta^* = (vartheta^*_1,..., vartheta^*_r)', with vartheta^*_i = sig^2_1 * rho_i + sig^2_2, have been estimated by REML, \hat{vartheta} = (\hat{vartheta}_1,...,\hat{vartheta}_r)', where \hat{vartheta}_i = U_i/nu_i.

Here, the goal is to test the following null hypotheses: H01: vartheta = vartheta^{H1} with vartheta^{H1}_i = sig^2_1 * rho_i + sig^2_2; H02: vartheta = vartheta^{H2} with vartheta^{H2}_i = 0 * rho_i + sig^2_2; H03: vartheta = vartheta^{H3} with vartheta^{H3}_i = 1 * rho_i + sig^2_2, for numerical values see Table 2.

clear 

## EXAMPLE - Unbalanced One-way random effects model

% Reset the global random number stream to its initial settings. s = RandStream('mt19937ar','Seed',0); RandStream.setGlobalStream(s); nGroups = 10; IncMatrix = 2 * (1:nGroups); model = OneWayModel(IncMatrix); 

## Set the TRUE model parameters

beta = 0; sigma1 = sqrt(.1); sigma2 = sqrt(1); 

## Generate the observations and/or the sufficient statistics

X = model.X; Z = model.Z; B = model.B; D = model.D; nObs = model.n; rho = model.rho; nu = model.nu; upsilon = sigma1 * randn(nGroups,1); epsilon = sigma2 * randn(nObs,1); nVar = numel(rho); y = X * beta + Z * upsilon + epsilon; t = B * y; u = zeros(numel(rho),1); for i = 1:numel(rho) u(i) = t' * D{i} * t; end varthetaREML = u./nu; 

## Compute the (REML) LRT test statistic

varthetaTRUE = sigma1^2*rho + sigma2^2; vartheta_H01 = varthetaTRUE; vartheta_H02 = 0*rho + sigma2^2; vartheta_H03 = 1*rho + sigma2^2; TABLE2 = [rho nu u varthetaTRUE varthetaREML vartheta_H01 vartheta_H02 vartheta_H03]; [~,id] = sort(rho,'descend'); TABLE2 = TABLE2(id,:); disp(TABLE2) lambda_1 = varthetaTRUE./vartheta_H01; lambda_2 = varthetaTRUE./vartheta_H02; lambda_3 = varthetaTRUE./vartheta_H03; LRT_H01 = (bsxfun(@minus,(u./varthetaTRUE).*lambda_1,nu)-bsxfun(@times,log(bsxfun(@rdivide,(u./varthetaTRUE).*lambda_1,nu)),nu)); LRT_H02 = (bsxfun(@minus,(u./varthetaTRUE).*lambda_2,nu)-bsxfun(@times,log(bsxfun(@rdivide,(u./varthetaTRUE).*lambda_2,nu)),nu)); LRT_H03 = (bsxfun(@minus,(u./varthetaTRUE).*lambda_3,nu)-bsxfun(@times,log(bsxfun(@rdivide,(u./varthetaTRUE).*lambda_3,nu)),nu)); LRT_H01 = sum(LRT_H01); LRT_H02 = sum(LRT_H02); LRT_H03 = sum(LRT_H03); disp([LRT_H01;LRT_H02;LRT_H03]) 
 Columns 1 through 7 19.2419 1.0000 0.6476 2.9242 0.6476 2.9242 1.0000 17.0439 1.0000 17.1167 2.7044 17.1167 2.7044 1.0000 14.8948 1.0000 2.7588 2.4895 2.7588 2.4895 1.0000 12.7658 1.0000 3.0109 2.2766 3.0109 2.2766 1.0000 10.6465 1.0000 0.4533 2.0646 0.4533 2.0646 1.0000 8.5313 1.0000 4.0207 1.8531 4.0207 1.8531 1.0000 6.4160 1.0000 0.5213 1.6416 0.5213 1.6416 1.0000 4.2959 1.0000 2.0575 1.4296 2.0575 1.4296 1.0000 2.1638 1.0000 0.8977 1.2164 0.8977 1.2164 1.0000 0.0000 100.0000 117.2464 1.0000 1.1725 1.0000 1.0000 Column 8 20.2419 18.0439 15.8948 13.7658 11.6465 9.5313 7.4160 5.2959 3.1638 1.0000 7.3095 18.7350 10.6475 

## Write the results into Table2.xls

xlswrite('Table2.xls',TABLE2); 

## Compute the LRT quantiles

prob = [.9;.95;.99]; q_LRT_H01 = lwChi2ConvInv(prob,lambda_1,nu,nu./lambda_1,nu./lambda_1,nu./lambda_1,'eynsham'); q_LRT_H02 = lwChi2ConvInv(prob,lambda_2,nu,nu./lambda_2,nu./lambda_2,nu./lambda_2,'eynsham'); q_LRT_H03 = lwChi2ConvInv(prob,lambda_3,nu,nu./lambda_3,nu./lambda_3,nu./lambda_3,'eynsham'); q_LRT_ASY = chi2inv(prob,nVar); disp([prob q_LRT_ASY q_LRT_H01 q_LRT_H02 q_LRT_H03]) 
 0.9000 15.9872 19.5778 25.5632 28.3141 0.9500 18.3070 22.2689 29.6080 31.2914 0.9900 23.2093 27.8874 38.6145 37.3990 

## Compute the exact LRT distributions

xmin = chi2inv(0.001,nVar); x_1 = linspace(xmin,q_LRT_H03(end),51)'; cdf_LRT_H01 = lwChi2ConvCdf(x_1,[],nu,nu.*(log(nu./lambda_1)-1),nu,lambda_1); x_2 = linspace(x_1(1),q_LRT_H03(end),51)'; cdf_LRT_H02 = lwChi2ConvCdf(x_2,[],nu,nu.*(log(nu./lambda_2)-1),nu,lambda_2); x_3 = linspace(x_1(1),q_LRT_H03(end),51)'; cdf_LRT_H03 = lwChi2ConvCdf(x_3,[],nu,nu.*(log(nu./lambda_3)-1),nu,lambda_3); cdf_LRT_ASY = chi2cdf(x_1,nVar); 

## Plot CDFs of the LRT test statistics for H01 H02 and H03

figure plot(x_1,cdf_LRT_ASY,'k:','LineWidth',1) grid hold on plot(x_1,cdf_LRT_H01,'b-','LineWidth',1) plot(x_2,cdf_LRT_H02,'m--','LineWidth',1) plot(x_3,cdf_LRT_H03,'r-.','LineWidth',1) hold off xlabel('lrt') ylabel('cdf') title('CDFs of the LRT statistics') 

## Create 'figure1.eps'

print -depsc figure1.eps