13. Design Matrices from Data and Model

Here, we will see how the design matrix can be constructed from the Data in a DataFrame and a model given as a string.

13.1. Data

using DataFrames, SparseArrays, LinearAlgebra
animal= ["animal1","animal2","animal3","animal4","animal5","animal6"]
sex   = ["m","f","f","m","f","f"]
breed = ["Angus","Angus","Hereford","Hereford","Angus","Angus"]
age   = [40,36,38,42,40,36]
df    = DataFrame(animal=animal,sex=sex,breed=breed,age=age,y=round.(randn(6),digits=3))

6 rows × 5 columns

	animal	sex	breed	age	y
	String	String	String	Int64	Float64
1	animal1	m	Angus	40	-0.473
2	animal2	f	Angus	36	-1.446
3	animal3	f	Hereford	38	-0.917
4	animal4	m	Hereford	42	1.196
5	animal5	f	Angus	40	-0.659
6	animal6	f	Angus	36	-0.294

13.2. Categorical variable

\[ y_{ij} = \mu + sex_i + e_{ij}. \]

We have seen previously how to construct the design matrix for a one-way model when the levels had sequential integer values. In this DataFrame the levels of sex are m and f. We will see below how to assign sequential integers to these strings “m” and “f”

13.2.1. Get the levels of sex from the DataFrame into a vector A

A=df[:,:sex]

6-element Vector{String}:
 "m"
 "f"
 "f"
 "m"
 "f"
 "f"

The unique function returns the unique levels of a vector:

res = unique(A)

2-element Vector{String}:
 "m"
 "f"

Now, we can make a dictionary where each unique level is the key and the associated value is the sequential integer:

dictA = Dict()                 # declare an empty dictionary
for (i,s) in enumerate(res)    # fill the dictionary with the values in res
    dictA[s] = i
end
dictA

Dict{Any, Any} with 2 entries:
  "f" => 2
  "m" => 1

Can use this dictionary to make design matrix:

ii = 1:size(A,1)            # get row numbers

jj = [dictA[i] for i in A]  # get column numbers using list comprehension

[ii jj]

6×2 Matrix{Int64}:
 1  1
 2  2
 3  2
 4  1
 5  2
 6  2

XA   = sparse(ii,jj,1.0)

6×2 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅ 
  ⋅   1.0
  ⋅   1.0
 1.0   ⋅ 
  ⋅   1.0
  ⋅   1.0

n = size(A,1)
Matrix([ones(n,1) XA])

6×3 Matrix{Float64}:
 1.0  1.0  0.0
 1.0  0.0  1.0
 1.0  0.0  1.0
 1.0  1.0  0.0
 1.0  0.0  1.0
 1.0  0.0  1.0

# This function returns a dictionary with the unique values in the vector "a" as the keys and their 
# sequential numbers as the associated values
# It also returns vector with the keys in sequential order. 
function mkDict(a)
  aUnique = unique(a)
  d = Dict()
  names = Array{String}(undef,size(aUnique,1))
  for (i,s) in enumerate(aUnique)
    names[i] = s
    d[s] = i
  end
  return d,names
end

mkDict (generic function with 1 method)

13.3. Continuous variable

\[ y_{ij} = \mu + sex_i + breed_j + \beta age_{ij} + e_{ij} \]

We will use the mkDict function to construct the design matrix for this model given the data in DataFrame `df’.

Get design matrix for sex

dictA,namesA   = mkDict(A)

(Dict{Any, Any}("f" => 2, "m" => 1), ["m", "f"])

namesA

2-element Vector{String}:
 "m"
 "f"

ii = 1:size(A,1)
jj = [dictA[i] for i in A]  #list comprehension 
[ii jj]
XA = sparse(ii,jj,1.0)

6×2 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅ 
  ⋅   1.0
  ⋅   1.0
 1.0   ⋅ 
  ⋅   1.0
  ⋅   1.0

Get design matrix for breed

B = df[:,:breed]
dictB,namesB   = mkDict(B)
jj   = [dictB[i] for i in B]  #list comprehension 
ii   = 1:size(B,1)
XB   = sparse(ii,jj,1.0)

6×2 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅ 
 1.0   ⋅ 
  ⋅   1.0
  ⋅   1.0
 1.0   ⋅ 
 1.0   ⋅ 

dictB

Dict{Any, Any} with 2 entries:
  "Angus"    => 1
  "Hereford" => 2

Get design matrix for age

CVal = df[:,:age]
CStr = fill("age",size(CVal,1))  # only one column in design matrix for age

6-element Vector{String}:
 "age"
 "age"
 "age"
 "age"
 "age"
 "age"

dictC,namesC   = mkDict(CStr)

(Dict{Any, Any}("age" => 1), ["age"])

jj   = [dictC[i] for i in CStr]  #list comprehension 
ii   = 1:size(CStr,1)

1:6

XC   = sparse(ii,jj,CVal)

6×1 SparseMatrixCSC{Int64, Int64} with 6 stored entries:
 40
 36
 38
 42
 40
 36

Get design matrix for all effects

n = size(A,1)
Matrix([ones(n,1) XA XB XC])

6×6 Matrix{Float64}:
 1.0  1.0  0.0  1.0  0.0  40.0
 1.0  0.0  1.0  1.0  0.0  36.0
 1.0  0.0  1.0  0.0  1.0  38.0
 1.0  1.0  0.0  0.0  1.0  42.0
 1.0  0.0  1.0  1.0  0.0  40.0
 1.0  0.0  1.0  1.0  0.0  36.0

["intercept"; namesA; namesB; namesC]

6-element Array{String,1}:
 "intercept"
 "m"        
 "f"        
 "Angus"    
 "Hereford" 
 "age"      

13.4. Two-way model with interaction

The \(\mathbf{X}\) matrix for the two-way model with interation between breed and sex

\[ y_{ijk} = \mu + sex_i + breed_j + sex_i*breed_j+ e_{ijk} \]

We already have the design matrices for the main effects.

13.4.1. Design Matrix for Interaction Term

13.4.1.1. Make vector of levels for interaction:

firstName = "Short"
lastName = " Course"
firstName * lastName

"Short Course"

AB = A .*" x ".*B

6-element Vector{String}:
 "m x Angus"
 "f x Angus"
 "f x Hereford"
 "m x Hereford"
 "f x Angus"
 "f x Angus"

Use AB to construct design matrix (XAB)

dictAB,namesAB   = mkDict(AB)
ii   = 1:size(A,1)
jj   = [dictAB[i] for i in AB]  #list comprehension 
XAB = sparse(ii,jj,1.0)

6×4 SparseMatrixCSC{Float64, Int64} with 6 stored entries:
 1.0   ⋅    ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅ 
  ⋅    ⋅   1.0   ⋅ 
  ⋅    ⋅    ⋅   1.0
  ⋅   1.0   ⋅    ⋅ 
  ⋅   1.0   ⋅    ⋅ 

namesAB

4-element Vector{String}:
 "m x Angus"
 "f x Angus"
 "f x Hereford"
 "m x Hereford"

13.4.1.2. Design Matrix

n = size(A,1)
Matrix([ones(n,1) XA XB XAB])

6×9 Matrix{Float64}:
 1.0  1.0  0.0  1.0  0.0  1.0  0.0  0.0  0.0
 1.0  0.0  1.0  1.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  0.0  1.0  0.0  0.0  1.0  0.0
 1.0  1.0  0.0  0.0  1.0  0.0  0.0  0.0  1.0
 1.0  0.0  1.0  1.0  0.0  0.0  1.0  0.0  0.0
 1.0  0.0  1.0  1.0  0.0  0.0  1.0  0.0  0.0

["intercept"; namesA; namesB; namesAB]

9-element Array{String,1}:
 "intercept"   
 "m"           
 "f"           
 "Angus"       
 "Hereford"    
 "m x Angus"   
 "f x Angus"   
 "f x Hereford"
 "m x Hereford"

13.5. Model with sex-specific slope for age

\[ y_{ij} = \mu + sex_i + \beta_{i}(age_{ij})+ e_{ij} \]

BVal = df[:,:age]
BStr = fill("age",size(BVal,1))

6-element Vector{String}:
 "age"
 "age"
 "age"
 "age"
 "age"
 "age"

AB = A.*" x ".*BStr 

6-element Vector{String}:
 "m x age"
 "f x age"
 "f x age"
 "m x age"
 "f x age"
 "f x age"

dAB,namesAB   = mkDict(AB)
ii    = 1:size(AB,1)
jj   = [dAB[i] for i in AB]  #list comprehension 
XAB   = sparse(ii,jj,BVal)
Matrix(XAB)

6×2 Matrix{Int64}:
 40   0
  0  36
  0  38
 42   0
  0  40
  0  36

13.5.1. Design Matrix for Model

Matrix([ones(n,1) XA XAB])

6×5 Matrix{Float64}:
 1.0  1.0  0.0  40.0   0.0
 1.0  0.0  1.0   0.0  36.0
 1.0  0.0  1.0   0.0  38.0
 1.0  1.0  0.0  42.0   0.0
 1.0  0.0  1.0   0.0  40.0
 1.0  0.0  1.0   0.0  36.0

["intercept"; namesA; namesAB]

5-element Vector{String}:
 "intercept"
 "m"
 "f"
 "m x age"
 "f x age"

13.6. Function to Construct Design Matrix

13.6.1. Function for main effects

We will do this by putting the code we have used earlier into a fucntion. First, let’s make our code work for quantitative or qualitiative factors.

# Test code for qualitative factor
factor = "sex"
cov    = false
data = df[:,Symbol(factor)]

n = size(data,1)
if cov==false
    str = data
    val = 1.0
else 
    str = fill(factor,n)
    val = data
end

dict,colNames   = mkDict(str)
ii    = 1:n                    # row numbers 
jj   = [dict[i] for i in str]  # column numbers
X    = sparse(ii,jj,val)
Matrix(X)    

6×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
 0.0  1.0
 1.0  0.0
 0.0  1.0
 0.0  1.0

# Test code for quantitative factor
factor = "age"
cov    = true

data = df[:,Symbol(factor)]
n = size(data,1)
if cov==false
    str = data
    val = 1.0
else 
    str = fill(factor,n)
    val = data
end

dict,colNames   = mkDict(str)
ii    = 1:n                    # row numbers 
jj   = [dict[i] for i in str]  # column numbers
X    = sparse(ii,jj,val)
Matrix(X)

6×1 Matrix{Int64}:
 40
 36
 38
 42
 40
 36

13.6.1.1. Put the code in a function:

function getX(factor,df;cov=false)
    data = df[:,Symbol(factor)]
    n = size(data,1)
    if cov==false
        str = data
        val = 1.0
    else 
        str = fill(factor,n)
        val = data
    end

    dict,names   = mkDict(str)
    ii    = 1:n                    # row numbers 
    jj   = [dict[i] for i in str]  # column numbers
    X    = sparse(ii,jj,val)    
end        

getX (generic function with 1 method)

X = getX("sex",df)
Matrix(X)

6×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
 0.0  1.0
 1.0  0.0
 0.0  1.0
 0.0  1.0

X = getX("age",df,cov=true)
Matrix(X)

6×1 Matrix{Int64}:
 40
 36
 38
 42
 40
 36

13.6.2. Function for main effects or interactions

function getX(factors,covs,df) #sex*age
    n = size(df,1)
    if covs[1] == false #sex
        str = df[:,Symbol(factors[1])]
        val = 1.0
    else
        str = fill(factors[1],n) 
        val = df[:,Symbol(factors[1])]    
    end       

    for i in 2:length(factors) #age
        if covs[i] == false
            str = str .*" x ".*df[:,Symbol(factors[i])]
            val = val .* 1.0 
        else
            str = str .*" x ".*fill(factors[i],n) 
            val = val .* df[:,Symbol(factors[i])]    
        end 
    end 
    dict,colNames   = mkDict(str)
    ii    = 1:n                    # row numbers 
    jj   = [dict[i] for i in str]  # column numbers
    X    = sparse(ii,jj,val)
    return X,colNames        
end            

getX (generic function with 2 methods)

Get design matrix for interaction between sex and breed

factors = ["sex", "breed"]
covs =[false, false]
X,colNames = getX(factors,covs,df)
Matrix(X)

6×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
 0.0  1.0  0.0  0.0
 0.0  1.0  0.0  0.0

colNames

4-element Vector{String}:
 "m x Angus"
 "f x Angus"
 "f x Hereford"
 "m x Hereford"

Get design matrix for interaction between age and age

factors = ["age", "age"]
covs =[true, true]
X,colNames = getX(factors,covs,df)
Matrix(X)

6×1 Matrix{Int64}:
 1600
 1296
 1444
 1764
 1600
 1296

colNames

1-element Array{String,1}:
 "age x age"

13.7. Get factors and covariables from model term

Consider model term: “sex*breed” and vector of covariables in the model: “age”

modelTerm   = "sex * breed"
covariables = ["age"];

split(modelTerm,"*")

2-element Vector{SubString{String}}:
 "sex "
 " breed"

factors = strip.(split(modelTerm,"*"))

2-element Vector{SubString{String}}:
 "sex"
 "breed"

covs = [i in covariables for i in factors]

2-element Vector{Bool}:
 0
 0

function getX(modelTerm,covariables,df)
    n = size(df,1)
    if modelTerm == "intercept"
        X = ones(n,1)
        colNames = ["intercept"]
        return X,colNames
    end
    factors = strip.(split(modelTerm,"*"))
    covs = [i in covariables for i in factors]
    
    if covs[1] == false
        str = string.(df[:,Symbol(factors[1])])
        val = 1.0
    else
        str = fill(factors[1],n) 
        val = df[:,Symbol(factors[1])]    
    end       

    for i in 2:length(factors)
        if covs[i] == false
            str = str .*" x ".*string.(df[:,Symbol(factors[i])])
            val = val .* 1.0 
        else
            str = str .*" x ".*fill(factors[i],n) 
            val = val .* df[:,Symbol(factors[i])]    
        end 
    end 
    dict,colNames   = mkDict(str)
    ii = 1:n                     # row numbers 
    jj = [dict[i] for i in str]  # column numbers
    X  = sparse(ii,jj,val)
    return X, strip(modelTerm)*": ".*colNames   
end

getX (generic function with 2 methods)

modelTerm   = "sex * breed"
covariables = ["age"]

X,colNames = getX(modelTerm,covariables,df)
Matrix(X)

6×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
 0.0  1.0  0.0  0.0
 0.0  1.0  0.0  0.0

colNames

4-element Vector{String}:
 "sex * breed: m x Angus"
 "sex * breed: f x Angus"
 "sex * breed: f x Hereford"
 "sex * breed: m x Hereford"

modelTerm = "breed*age"
X,colNames = getX(modelTerm,covariables,df)
Matrix(X)

6×2 Matrix{Float64}:
 40.0   0.0
 36.0   0.0
  0.0  38.0
  0.0  42.0
 40.0   0.0
 36.0   0.0

colNames

2-element Vector{String}:
 "breed*age: Angus x age"
 "breed*age: Hereford x age"

modelTerm = "sex"
X,colNames = getX(modelTerm,covariables,df)
Matrix(X)

6×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0
 0.0  1.0
 1.0  0.0
 0.0  1.0
 0.0  1.0

colNames

2-element Vector{String}:
 "sex: m"
 "sex: f"

modelTerm = "intercept"
X,colNames = getX(modelTerm,covariables,df)
Matrix(X)

6×1 Matrix{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

colNames

1-element Vector{String}:
 "intercept"

13.8. Get model-terms from string representaion of model

Consider model:

modelEq = "y = intercept + sex + breed + sex*breed + age"

"y = intercept + sex + breed + sex*breed + age"

modelParts = strip.(split(modelEq,"="))

2-element Vector{SubString{String}}:
 "y"
 "intercept + sex + breed + sex*breed + age"

depVar = modelParts[1]
modelTerms = strip.(split(modelParts[2],"+"))

5-element Vector{SubString{String}}:
 "intercept"
 "sex"
 "breed"
 "sex*breed"
 "age"

X,colNames = getX(modelTerms[1],covariables,df)
Matrix(X)

6×1 Matrix{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0
 1.0

i = 2
Xi,namesi = getX(modelTerms[i],covariables,df)
X = [X Xi]
colNames = [colNames; namesi]
Matrix(X)

6×3 Matrix{Float64}:
 1.0  1.0  0.0
 1.0  0.0  1.0
 1.0  0.0  1.0
 1.0  1.0  0.0
 1.0  0.0  1.0
 1.0  0.0  1.0

colNames

3-element Vector{String}:
 "intercept"
 "sex: m"
 "sex: f"

i = 3
Xi,namesi = getX(modelTerms[i],covariables,df)
X = [X Xi]
colNames = [colNames; namesi]
Matrix(X)

6×5 Matrix{Float64}:
 1.0  1.0  0.0  1.0  0.0
 1.0  0.0  1.0  1.0  0.0
 1.0  0.0  1.0  0.0  1.0
 1.0  1.0  0.0  0.0  1.0
 1.0  0.0  1.0  1.0  0.0
 1.0  0.0  1.0  1.0  0.0

colNames

5-element Vector{String}:
 "intercept"
 "sex: m"
 "sex: f"
 "breed: Angus"
 "breed: Hereford"

function getLhsRhs(modelEq,covariables,df)
    modelParts = strip.(split(modelEq,"="))
    depVar = modelParts[1]
    y = df[:,Symbol(depVar)]
    modelTerms = strip.(split(modelParts[2],"+"))
    X,colNames = getX(modelTerms[1],covariables,df)
    for i = 2:size(modelTerms,1)
        Xi,namesi = getX(modelTerms[i],covariables,df)
        X = [X Xi]
        colNames = [colNames; namesi]
    end
    return X'X,X'y,colNames
end

getLhsRhs (generic function with 1 method)

lhs,rhs,colNames,X = getLhsRhs(modelEq,covariables,df)
[colNames Matrix(lhs) rhs]

10×12 Matrix{Any}:
 "intercept"                  6.0  …    3.0   1.0   1.0   232.0   -2.593
 "sex: m"                     2.0       0.0   0.0   1.0    82.0    0.723
 "sex: f"                     4.0       3.0   1.0   0.0   150.0   -3.316
 "breed: Angus"               4.0       3.0   0.0   0.0   152.0   -2.872
 "breed: Hereford"            2.0       0.0   1.0   1.0    80.0    0.279
 "sex*breed: m x Angus"       1.0  …    0.0   0.0   0.0    40.0   -0.473
 "sex*breed: f x Angus"       3.0       3.0   0.0   0.0   112.0   -2.399
 "sex*breed: f x Hereford"    1.0       0.0   1.0   0.0    38.0   -0.917
 "sex*breed: m x Hereford"    1.0       0.0   0.0   1.0    42.0    1.196
 "age: age"                 232.0     112.0  38.0  42.0  9000.0  -92.534

Matrix(lhs)

10×10 Matrix{Float64}:
   6.0   2.0    4.0    4.0   2.0   1.0    3.0   1.0   1.0   232.0
   2.0   2.0    0.0    1.0   1.0   1.0    0.0   0.0   1.0    82.0
   4.0   0.0    4.0    3.0   1.0   0.0    3.0   1.0   0.0   150.0
   4.0   1.0    3.0    4.0   0.0   1.0    3.0   0.0   0.0   152.0
   2.0   1.0    1.0    0.0   2.0   0.0    0.0   1.0   1.0    80.0
   1.0   1.0    0.0    1.0   0.0   1.0    0.0   0.0   0.0    40.0
   3.0   0.0    3.0    3.0   0.0   0.0    3.0   0.0   0.0   112.0
   1.0   0.0    1.0    0.0   1.0   0.0    0.0   1.0   0.0    38.0
   1.0   1.0    0.0    0.0   1.0   0.0    0.0   0.0   1.0    42.0
 232.0  82.0  150.0  152.0  80.0  40.0  112.0  38.0  42.0  9000.0

rhs

10-element Vector{Float64}:
  -2.5930000000000004
   0.723
  -3.3160000000000003
  -2.8720000000000003
   0.2789999999999999
  -0.473
  -2.399
  -0.917
   1.196
 -92.534

QRLhs = qr(lhs) 
sol = QRLhs\rhs

10-element Vector{Float64}:
 -2.921499999996586
  1.902000000000359
  0.0
  0.1524999999999361
  0.0
 -1.7160000000001157
  0.0
  0.0
  0.0
  0.05274999999991032