MATLAB Examples

Contents

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