Question regarding implementation of Multivariate Generalized Linear Model
2 views (last 30 days)
Show older comments
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
Answers (0)
See Also
Categories
Find more on Copula Distributions and Correlated Samples in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!