## Still can't make it work [🇷 for BE/BA]

❝ Here all subjects are (y is logarithmised)

# Theta estimate θ was taken by REML

using CSV, DataFrames, StatsModels, StatsBase path = dirname(@__FILE__) cd(path) data   = CSV.File("dat.csv", delim=' ')|> DataFrame #Sort data sort!(data, [:Subj, :Trt, :Per]) #Get X matrix, Z, y X = ModelMatrix(ModelFrame(@formula(0 ~ Seq + Per + Trt), data)).m Z = ModelMatrix(ModelFrame(@formula(0 ~ 0 + Trt), data, contrasts = Dict(:Trt => StatsModels.FullDummyCoding()))).m y   = data[!, :y] #Xv vector of Xi, Zv vector of Zi, yv vector of yi u = unique(data[!, :Subj]) Xv = Vector{Matrix}(undef, length(u)) Zv = Vector{Matrix}(undef, length(u)) yv = Vector{Vector}(undef, length(u)) for i = 1:length(u)     v = findall(x -> x == u[i], data[!, :Subj])     Xv[i] = view(X, v, :)     Zv[i] = view(Z, v, :)     yv[i] = view(y, v) end # Theta estimate θ[1:2] for R, θ[1:3] for G # Very hard to take good θ estmate for this design # If you provide your estimate, β can be recalculated θ = [0.013246492714940418, 0.008891008058562478, 0.03621599611178057, 0.06160355780666301, 0.9661995154179528] #G matrix G = [θ[3] sqrt(θ[3]*θ[4])*θ[5]; sqrt(θ[3]*θ[4])*θ[5] θ[4]] #Vector of R matrices Rv = Diagonal.(map(x -> x * θ[1:2], Zv)) #Construct vector of Vi Vv = Vector{Matrix}(undef, length(u)) for i = 1:length(u)     global Vv[i] = Zv[i]*G*Zv[i]' + Rv[i] end #Vector of inverted Vi iVv = inv.(Vv) M1 = zeros(6, 6) M2 = zeros(6) #Calc M1 & M2 for i = 1:length(u)     global M1 .+= Xv[i]'*iVv[i]*Xv[i]     global M2 .+= Xv[i]'*iVv[i]*yv[i] end β = inv(M1) * M2 #= julia> β 6-element Array{Float64,1}:  7.904258681084915  0.0547761264037151  0.05092362547466389  0.0012959346740553102  0.048118895192829976  0.02239133333333365 =#

I get beta:

7.904258681084915
0.0547761264037151
0.05092362547466389
0.0012959346740553102
0.048118895192829976
0.02239133333333365