Ntrials=10; Nsub=12; N=Ntrials*Nsub; pop_effect=10; sigma_b_2=1; ymin=pop_effect-3*sqrt(sigma_b_2); ymax=pop_effect+3*sqrt(sigma_b_2); % Within subject variance lambda=rand(Nsub,1); %lambda(5)=10; % Generate design matrices X1=kron(eye(Nsub),ones(Ntrials,1)); X2=ones(Nsub,1); % Generate first-level covariance components C1=zeros(N,N); for n=1:Nsub, q=zeros(Ntrials*Nsub,1); start=(n-1)*Ntrials+1; stop=start+Ntrials-1; q(start:stop)=1; Q{n}=diag(q); C1=C1+lambda(n)*Q{n}; end % Second level lambda(n+1)=sigma_b_2; C2=lambda(n+1)*eye(Nsub); Q{Nsub+1}=X1*eye(Nsub)*X1'; % Error covariance for collapsed model C=C1+X1*C2*X1'; % Generate and plot data sub=randn(Nsub,1)*sqrt(lambda(Nsub+1))+pop_effect; y=[]; figure for n=1:Nsub, trials=randn(Ntrials,1)*sqrt(lambda(n))+sub(n); subplot(1,Nsub,n); plot(ones(Ntrials,1),trials,'x'); axis([0 2 ymin ymax]); set(gca,'XTickLabel',[]); if n>1 set(gca,'YTickLabel',[]); end xlabel(int2str(n)); y=[y;trials]; end % Plot design matrix and covariance components figure subplot(2,2,1); imagesc(X1); colormap gray title('X_1'); subplot(2,2,2); imagesc(C1); colormap gray title('C_1'); subplot(2,2,3); imagesc(X2); colormap gray title('X_2'); subplot(2,2,4); imagesc(C2); colormap gray title('C_2'); % ReML estimation X=X1*X2; YY=y*y'; [RC,Rh,RPh]=spm_reml(YY,X,Q,1); % Display results figure true_var=lambda(1:Nsub); est_var=Rh(1:Nsub); mv=max([true_var; est_var]); plot(true_var,est_var,'rx'); hold on plot([0 mv],[0 mv],'-'); title('Error variances'); xlabel('True'); ylabel('ReML estimate');