spm_ancova.m 2.9 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889
  1. function [F,df,beta,xX,xCon] = spm_ancova(xX,V,Y,c)
  2. % estimation and inference of a linear model
  3. % FORMAT [F,df,beta,xX,xCon] = spm_ancova(xX,V,Y,c);
  4. %
  5. % xX - (m x p) Design matrix or structure
  6. % V - (m x m) error covariance constraint
  7. % Y - {m x n} matrix of response {m x 1} variables
  8. % c - {p x q} matrix of (q) contrasts
  9. %
  10. % F - {t x n} matrix of T or F values
  11. % df - {1 x 2} vector of degrees of freedom
  12. % beta - {p x n} matrix of parameter estimates
  13. % xX - design matrix structure
  14. % xCon - contrast structure
  15. %__________________________________________________________________________
  16. %
  17. % spm_ancova uses a General Linear Model of the form:
  18. %
  19. % Y = X*beta + K*e
  20. %
  21. % to compute the parameter estimates (beta) and make inferences (T or F)
  22. % where V = K*K' represents the correlation structure. If c has only one
  23. % column T statistics are returned, otherwise F ratios are computed.
  24. %__________________________________________________________________________
  25. % Copyright (C) 2002-2011 Wellcome Trust Centre for Neuroimaging
  26. % Karl Friston
  27. % $Id: spm_ancova.m 7033 2017-03-05 11:19:18Z karl $
  28. % assume V = I (i.i.d.) if empty
  29. %--------------------------------------------------------------------------
  30. if isempty(V)
  31. V = speye(size(Y,1));
  32. end
  33. % or scalar
  34. %--------------------------------------------------------------------------
  35. if isscalar(V)
  36. V = speye(size(Y,1))*V;
  37. end
  38. % create design matrix structure if necessary
  39. %--------------------------------------------------------------------------
  40. if ~isstruct(xX)
  41. xX = spm_sp('Set',xX);
  42. end
  43. if ~isfield(xX,'pX')
  44. xX.pX = spm_sp('x-',xX);
  45. end
  46. % estimate parameters and sum of squared residuals
  47. %--------------------------------------------------------------------------
  48. beta = xX.pX*Y; %-Parameter estimates
  49. res = spm_sp('r',xX,Y); %-Residuals
  50. ResSS = sum(res.^2); %-Res sum-of-squares
  51. % default contrast
  52. %--------------------------------------------------------------------------
  53. if nargin < 4
  54. c = speye(size(beta,1));
  55. end
  56. % contrast
  57. %--------------------------------------------------------------------------
  58. xCon = spm_FcUtil('Set','','F','c',c,xX);
  59. h = spm_FcUtil('Hsqr',xCon,xX);
  60. X1o = spm_FcUtil('X1o',xCon,xX);
  61. % ensure trace(V) = m and get traces
  62. %--------------------------------------------------------------------------
  63. V = V*length(V)/trace(V);
  64. [trRV, trRVRV] = spm_SpUtil('trRV',xX,V);
  65. [trMV, trMVMV] = spm_SpUtil('trMV',X1o,V);
  66. df = [trMV^2/trMVMV trRV^2/trRVRV];
  67. % F statistics
  68. %--------------------------------------------------------------------------
  69. F = sum((h*beta).^2,1)./(ResSS*trMV/trRV);
  70. if size(c,2) == 1
  71. % T statistics
  72. %----------------------------------------------------------------------
  73. F = sqrt(F).*sign(c'*beta);
  74. end