Makes a fixed and random effects term into a mixed effects model.
Usage: model = random( ran [,fix [,strran [,strfix [,ranisvar]]] );
Internally a random consists of two terms that describe linear models for
the mean and the variance. The variables for the variance term are
symmetric matrices stored as column vectors. If n is the number of rows
in the mean term, then n^2 is the number of rows in the variance term.
ran = term, or anything that can be converted into a term, for the
variance. The single variable for the variance term is
vec(double(ran)*double(ran)'). If ran is a scalar s, then it is
vec(identity matrix)*s^2 whose size matches that of the other
term in a binary operator, and its name is 'Is^2' or 'I' if s=1.
fix = term, or anything that can be converted into a term, for the
mean, empty if empty or absent. ran and fix must have the same
number of rows, except when one is a scalar.
strran = string for name of ran. If empty or absent, 'ran' is used. If
ran is a scalar, then 'I' is used.
strfix = cell array of strings for names for the variables in fix. If
empty or absent, chosen by term(fix).
ranisvar = 1 if ran is already a term for the variance.
If ran is a random, model=ran. With no arguments, model is empty.
Let model.mean be the mean term, and model.variance be the variance term.
The following operators are overloaded for random models m1 and m2:
m = m1 + m2: m.mean = m1.mean + m2.mean, and
m.variance = m1.variance + m2.variance.
m = m1 - m2: m.mean = m1.mean - m2.mean, and
m.variance = m1.variance - m2.variance.
m = m1 * m2: m.mean = m1.mean * m2.mean, and
m.variance = m1.variance * m2.variance +
m1.variance * sum_{i<=j} random( (m2.mean_i+m2.mean_j)/2 ) +
m1.variance * sum_{i<j} random( (m2.mean_i-m2.mean_j)/2 ) +
sum_{i<=j} random( (m1.mean_i+m1.mean_j)/2 ) * m2.variance +
sum_{i<j} random( (m1.mean_i-m1.mean_j)/2 ) * m2.variance,
where _i indicates variable i. m_.mean_i is divided by its
maximum absolute value so the variables have comparable
size, and the factor of 1/2 produces nicer printout.
These last four extra sums are equivalent to
sum_{i<=j} vec( m1.mean_i * m1.mean_j' ) * m2.variance +
sum_{i<=j} vec( m2.mean_i * m2.mean_j' ) * m1.variance,
which allow for an interaction between fixed effects and
random effects with arbitrary variance matrices that is
invariant to any linear transformation of the fixed effects.
The advantage of the former parameterisation is that if its
coefficients are positive then the variance matrix of the
observations is positive definite. The reason is that
element-wise products of positive definite matrices are
positive definite. However enforcing positivity has two bad
side effects: we lose invariance to any linear
transformation of the fixed effects, and it restricts the
range of positive definite matrices in the model.
Specifically, if r is the correlation between two fixed
effects (crossed with a random effect) with standard
deviations sd1<=sd2, then |r|<=sd1/sd2. In other words, if
the sd's are very different, then the correlaton can't be
too big.
m = m1 ^ k: product of m1 with itself, k times.
Algebra: commutative, associative and distributive rules for random a,b,c
a + b = b + a
a * b = b * a
(a + b) + c = a + (b + c)
(a * b) * c = a * (b * c)
(a + b) * c = a * c + b * c only if a,b,c are purely fixed or random
a + 0 = a
a * 1 = a
a * 0 = 0
However note that
a + b - c =/= a + (b - c) =/= a - c + b
If t1,t2 are terms:
random(t1 * t2) = random(t1) * random(t2), but
random(t1 + t2) =/= random(t1) + random(t2), since the LHS is one term
The following functions are overloaded for random model:
char(model) = [char(model.mean), char(model.variance)].
double(model) = [double(model.mean), double(model.variance)].
size(model) = [size(model.mean), size(model.variance)].
isempty(model) = 1 if model is empty and 0 if not.0001 function model = random( ran, fix, strran, strfix, ranisvar ); 0002 0003 %Makes a fixed and random effects term into a mixed effects model. 0004 % 0005 % Usage: model = random( ran [,fix [,strran [,strfix [,ranisvar]]] ); 0006 % 0007 % Internally a random consists of two terms that describe linear models for 0008 % the mean and the variance. The variables for the variance term are 0009 % symmetric matrices stored as column vectors. If n is the number of rows 0010 % in the mean term, then n^2 is the number of rows in the variance term. 0011 % 0012 % ran = term, or anything that can be converted into a term, for the 0013 % variance. The single variable for the variance term is 0014 % vec(double(ran)*double(ran)'). If ran is a scalar s, then it is 0015 % vec(identity matrix)*s^2 whose size matches that of the other 0016 % term in a binary operator, and its name is 'Is^2' or 'I' if s=1. 0017 % fix = term, or anything that can be converted into a term, for the 0018 % mean, empty if empty or absent. ran and fix must have the same 0019 % number of rows, except when one is a scalar. 0020 % strran = string for name of ran. If empty or absent, 'ran' is used. If 0021 % ran is a scalar, then 'I' is used. 0022 % strfix = cell array of strings for names for the variables in fix. If 0023 % empty or absent, chosen by term(fix). 0024 % ranisvar = 1 if ran is already a term for the variance. 0025 % If ran is a random, model=ran. With no arguments, model is empty. 0026 % 0027 % Let model.mean be the mean term, and model.variance be the variance term. 0028 % The following operators are overloaded for random models m1 and m2: 0029 % 0030 % m = m1 + m2: m.mean = m1.mean + m2.mean, and 0031 % m.variance = m1.variance + m2.variance. 0032 % 0033 % m = m1 - m2: m.mean = m1.mean - m2.mean, and 0034 % m.variance = m1.variance - m2.variance. 0035 % 0036 % m = m1 * m2: m.mean = m1.mean * m2.mean, and 0037 % m.variance = m1.variance * m2.variance + 0038 % m1.variance * sum_{i<=j} random( (m2.mean_i+m2.mean_j)/2 ) + 0039 % m1.variance * sum_{i<j} random( (m2.mean_i-m2.mean_j)/2 ) + 0040 % sum_{i<=j} random( (m1.mean_i+m1.mean_j)/2 ) * m2.variance + 0041 % sum_{i<j} random( (m1.mean_i-m1.mean_j)/2 ) * m2.variance, 0042 % 0043 % where _i indicates variable i. m_.mean_i is divided by its 0044 % maximum absolute value so the variables have comparable 0045 % size, and the factor of 1/2 produces nicer printout. 0046 % These last four extra sums are equivalent to 0047 % 0048 % sum_{i<=j} vec( m1.mean_i * m1.mean_j' ) * m2.variance + 0049 % sum_{i<=j} vec( m2.mean_i * m2.mean_j' ) * m1.variance, 0050 % 0051 % which allow for an interaction between fixed effects and 0052 % random effects with arbitrary variance matrices that is 0053 % invariant to any linear transformation of the fixed effects. 0054 % The advantage of the former parameterisation is that if its 0055 % coefficients are positive then the variance matrix of the 0056 % observations is positive definite. The reason is that 0057 % element-wise products of positive definite matrices are 0058 % positive definite. However enforcing positivity has two bad 0059 % side effects: we lose invariance to any linear 0060 % transformation of the fixed effects, and it restricts the 0061 % range of positive definite matrices in the model. 0062 % Specifically, if r is the correlation between two fixed 0063 % effects (crossed with a random effect) with standard 0064 % deviations sd1<=sd2, then |r|<=sd1/sd2. In other words, if 0065 % the sd's are very different, then the correlaton can't be 0066 % too big. 0067 % 0068 % m = m1 ^ k: product of m1 with itself, k times. 0069 % 0070 % Algebra: commutative, associative and distributive rules for random a,b,c 0071 % a + b = b + a 0072 % a * b = b * a 0073 % (a + b) + c = a + (b + c) 0074 % (a * b) * c = a * (b * c) 0075 % (a + b) * c = a * c + b * c only if a,b,c are purely fixed or random 0076 % a + 0 = a 0077 % a * 1 = a 0078 % a * 0 = 0 0079 % However note that 0080 % a + b - c =/= a + (b - c) =/= a - c + b 0081 % If t1,t2 are terms: 0082 % random(t1 * t2) = random(t1) * random(t2), but 0083 % random(t1 + t2) =/= random(t1) + random(t2), since the LHS is one term 0084 % 0085 % The following functions are overloaded for random model: 0086 % char(model) = [char(model.mean), char(model.variance)]. 0087 % double(model) = [double(model.mean), double(model.variance)]. 0088 % size(model) = [size(model.mean), size(model.variance)]. 0089 % isempty(model) = 1 if model is empty and 0 if not. 0090 0091 if nargin==0 0092 model.mean=[]; 0093 model.variance=[]; 0094 model=class(model,'random'); 0095 elseif isa(ran,'random') 0096 model=ran; 0097 else 0098 model.mean=term; 0099 model.variance=term; 0100 if ~isempty(ran) 0101 if nargin<3 || isempty(strran) 0102 strran=inputname(1); 0103 end 0104 if ~isa(ran,'term') 0105 ran=term(ran,strran); 0106 end 0107 if nargin==5 & ranisvar 0108 model.variance=ran; 0109 else 0110 x=double(ran); 0111 if prod(size(ran))==1 0112 if x==1 0113 strran='I'; 0114 else 0115 strran=['I' num2str(x) '^2']; 0116 end 0117 end 0118 v=x*x'; 0119 model.variance=term(v(:),strran); 0120 end 0121 end 0122 if nargin>=2 && ~isempty(fix) 0123 if nargin<4 || isempty(strfix) 0124 strfix=inputname(2); 0125 end 0126 if ~isa(fix,'term') 0127 fix=term(fix,strfix); 0128 end 0129 model.mean=fix; 0130 end 0131 model=class(model,'random'); 0132 end 0133 superiorto('term'); 0134 0135 return 0136 end