how to compute a linear mixed effect using nlmefit?

Question

Cyril Pernet 2013-3-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67041-how-to-compute-a-linear-mixed-effect-using-nlmefit

Hi,

I'd like to fit a simple linear mixed model Y=XB+Zb+e with X the design matrix of the fixed effect (always the same size but not always the same values) and Z described subjects (which I specify as random)

here is the code I used

% generate data
for subject=1:10
   x(:,subject) = [1:10]+randi(30,1);
   coef(subject) = (3+rand(1));
   y(:,subject) = coef(subject)*x(:,subject)+3*randn(10,1)- 5*mean(x(:,subject));
end
% create X, Y, subject for nlmefit
Y = y(:);
X = [x(:) ones(100,1)]; 
subject = sum(kron(diag(1:10),ones(10,1)),2);
% fit the data using 'model'
model = @(Betas,X) (X*Betas)
[Betas,PSI,stats] = nlmefit(X,Y,subject,[],model,[1 0])

The error is: Error using * Inner matrix dimensions must agree. Error in @(Betas,X)(X*Betas)

in a fixed effect Betas=pinv(X)Y and the fitted data = X*Betas, and that why i defined model this way, assuming that for each subject, parameters are fitted using 'model' ?? any idea what I am doing wrong ?

Thanks Cyril

2 个评论
显示无隐藏无

Walter Roberson 2013-3-13

Do you in fact want algebraic matrix multiplication? Or do you want element-by-element multiplication which is the .* operator ?

Cyril Pernet 2013-3-14

在 MATLAB Online 中打开

I'm looking for matrix multiplication (for fixed effect Betas = pinv(X)*Y and thus model = X*Betas) - I thought maybe using twice the same variable name confused matlab so I changed the code to

model = @(betas,X) (X*betas);
beta0 = [1 0]';
[Betas,PSI,stats] = nlmefit(X,Y,subject,[],model,beta0)

since the error was dimension issue - now for sure X*beta0 (the 1st guess for Betas) works - still I have the same error message :-(

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Tom Lane 2013-3-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/67041-how-to-compute-a-linear-mixed-effect-using-nlmefit#answer_78699

在 MATLAB Online 中打开

It looks like nlmefit invokes your model function with betas as a row vector. Try this:

model = @(Betas,X) (X*Betas(:))

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Cyril Pernet 2013-3-15

编辑：Cyril Pernet 2013-3-15

在 MATLAB Online 中打开

oops, this was indeed the issue - thx

now I got this, I'd like to generate some F test, I used the following code

Res = reshape(stats.ires,10,10);
for s=1:10
  Yhat(:,s) = y(:,s) - Res(:,s);
end
SSeffect = norm(Yhat(:)-mean(Yhat(:))).^2;
SStotal = norm(Y-mean(Y)).^2;
R2 = SSeffect / SStotal;
SSerror = norm(Res(:)-mean(Res(:))).^2;
df = (rank(X)-1)+9; % add number of subjects -1? 
                  % alternatively it can be nb of subjects - rank(X)
dfe = stats.dfe; 
F = (SSeffect/df) / (SSerror/dfe);
p_val = 1-fcdf(F,df,dfe);

obviously that works, doesn't mean that's right! in particular I'm not sure

(1) if the modeled data correspond to the data minus the individual residuals (stats.ires) or the weighted residuals (stats.cwres)

(2) of the df of the model ? for each subject we fit X so df is rank(X)-1, we also fit each subject, so maybe add the number of subjects -1 ? - alternative;y, I read the is can be the number of subjects - rank(X)

No worries, if you don't know this, that's already a relieve to get the code running ..

请先登录，再进行评论。