Asked by 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

Related Content

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

Learn moreOpportunities for recent engineering grads.

Apply Today
## 0 Comments