Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

New to MATLAB?

Question regarding implementation of Multivariate Generalized Linear Model

Asked by Won Hwa Kim

Won Hwa Kim

on 2 Apr 2013

I am trying to implement Multivariate General Linear Model in Matlab, which is implemented in R.

When I put in some random input in both R and my implementation, I get different answers and I am not sure how to proceed.

I know where I am doing it wrong, but I cannot fix it since I do not know the details of the theory.

It would be really great if anyone can tell me what I am doing wrong here.

Thanks.

if true
% MGLM using Wilk's Lambda
% Model: Y = XB + E
%
% Input
%     Y: response multivariate variable, n x m
%     X: model matrix, n x (k+1)
%     
% Output
%     stats.F: F-stats
%     stats.p: p-value
%
clc; clear;
%% Variables
mu1= [0 1 2; 0 2 2; 0 3 2;0 4 3]; % group 1
mu2=[2 0 1; 0 0 1;1 0 1;3 0 1]; % group 2
Y = [mu1;mu2];  % response variable
X = [1 1 1 1 0 0 0 0]'; % model 
%% dimensions
n=size(Y,1);
k=size(X,2); % Constant
p=size(X,2);
%% Computing sum-of-squares
temp = inv([ones(n,1) X]'*[ones(n,1) X])*[ones(n,1) X]'*Y;
intercept = temp(1,:);  % intercept
eB = temp(2,:); % B estimation
eY = X*eB;  % Y estimation
eE = Y - repmat(intercept,n,1) - eY;    % residual
Y2 = Y - repmat(intercept,n,1);
SSP_R = eE'*eE; % SSP residual matrix
y_mean = mean(Y)';  % mean of y
SSP_Reg = eY' * eY - n*y_mean*y_mean';
% Where I am doing it wrong....
L = ones(n,1);  % Hypothesis matrix L 
SSP_H = (eB'*L')*pinv(L*pinv(X'*X)*L')*(L*eB); 
%% Wilk's lambda 
WL = det(SSP_R) / det(SSP_H + SSP_R);
%% F-stat and p-value approximation 
m = size(SSP_R,1); 
s = p; 
r = n-k-1 - (m-s+1)/2;
u = (m*s-2)/4;
if m^2 + s^2 - 5 > 0
  t = sqrt((m^2 * s^2 - 4) / (m^2 + s^2 - 5));
else
  t = 0;
end
F = (1-WL^(1/t))/ WL^(1/t) * (r*t - 2*u)/(m*s);
pval = 1-fcdf(F,m*s, r*t-2*u);
%% Output
stats.F = F;
stats.p = pval;
display(stats);
end

0 Comments

Won Hwa Kim

Won Hwa Kim

Products

0 Answers

Contact us