Construct Design Matrix for Regression on Matrix I have a model 
$Y_{i,j} = a_i + b_j$,
where Y is a NxK matrix, i.e. a is a vector of size N and b is a vector of size K. If i want to write this model as a linear regression equation, i.e. with a Design matrix $X$ ($Y=X\beta + e)$, I seem to get into trouble. My design matrix looks like a 3 dimensional tensor and I have to introduce contraction over the third axis to get my formulation to work. I think I am doing something wrong here. How would I need to proceed to construct the design matrix in this case?
 A: Here is how you can do an example using R. First making factors for rows and columns, and then calling model.matrix:
     rows <- rep(1:3, 3)
     rows
    [1] 1 2 3 1 2 3 1 2 3
     cols <- rep(1:3, rep(3, 3))
     cols
    [1] 1 1 1 2 2 2 3 3 3
     rows <- as.factor(rows) ; cols <- as.factor(cols)

Then making the design matrix, omitting the intercept (which we can do, since intercept is sum of all rows, also is sum of all columns):
     X <- model.matrix(~ 0+rows+cols)
     X
      rows1 rows2 rows3 cols2 cols3
    1     1     0     0     0     0
    2     0     1     0     0     0
    3     0     0     1     0     0
    4     1     0     0     1     0
    5     0     1     0     1     0
    6     0     0     1     1     0
    7     1     0     0     0     1
    8     0     1     0     0     1
    9     0     0     1     0     1
    attr(,"assign")
    [1] 1 1 1 2 2
    attr(,"contrasts")
    attr(,"contrasts")$rows
    [1] "contr.treatment"
    
    attr(,"contrasts")$cols
    [1] "contr.treatment"

