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
Complete thread:
- Semireplicated + REML in R ElMaestro 2020-07-10 13:19 [🇷 for BE/BA]
- Semireplicated + REML in R Helmut 2020-07-10 19:03
- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57
- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43
- Braveheart! Helmut 2020-07-11 15:38
- Who can help? ElMaestro 2020-07-12 12:28
- Who can help? ElMaestro 2020-07-12 13:26
- Update II ElMaestro 2020-07-12 21:36
- Update III ElMaestro 2020-07-12 21:46
- Final update today ElMaestro 2020-07-12 22:27
- Medium rare. Helmut 2020-07-13 13:52
- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18
- Will take much more hours still… Helmut 2020-07-13 15:34
- Negative determinant ElMaestro 2020-07-14 03:22
- Are we loosers? Helmut 2020-07-14 13:58
- "we"? Loosers? ElMaestro 2020-07-14 15:07
- Misunderstanding? Helmut 2020-07-14 15:32
- "we"? Loosers? ElMaestro 2020-07-14 15:07
- Are we loosers? Helmut 2020-07-14 13:58
- took just 52 hrs to do it :-) ElMaestro 2020-07-13 14:18
- Medium rare. Helmut 2020-07-13 13:52
- Braveheart! ElMaestro 2020-07-13 10:13
- Braveheart! Helmut 2020-07-13 14:16
- Braveheart! PharmCat 2020-07-15 14:19
- Braveheart! ElMaestro 2020-08-02 17:39
- Who can help? ElMaestro 2020-07-12 12:28
- Braveheart! Helmut 2020-07-11 15:38
- Avoid partial replicate designs, pleeeze! ElMaestro 2020-07-11 14:43
- Avoid partial replicate designs, pleeeze! Helmut 2020-07-11 11:57
- Semireplicated + REML in R ElMaestro 2020-07-10 19:24
- We were all blind (except Detlew) Helmut 2020-07-15 14:27
- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- Desultory thoughts Helmut 2020-07-15 17:33
- FDA RSABE is ISC d_labes 2020-07-15 18:13
- FDA RSABE is ISC Helmut 2020-07-16 11:11
- Desultory thoughts ElMaestro 2020-07-15 23:06
- Desultory thoughts Helmut 2020-07-16 10:59
- FDA RSABE is ISC d_labes 2020-07-15 18:13
- Desultory thoughts Helmut 2020-07-15 17:33
- Phoenix - which template? mittyri 2020-07-19 00:42
- FDA RSABE Project template_ v1.4.phxproj Helmut 2020-07-19 01:45
- It is the opposite way around for me ElMaestro 2020-07-15 16:25
- "By popular demand": likelihood ElMaestro 2020-07-24 10:07
- And by the way.... ElMaestro 2020-07-24 12:52
- And by the way.... PharmCat 2020-08-03 14:24
- Not understood ElMaestro 2020-08-03 22:55
- Not understood PharmCat 2020-08-05 01:41
- Not understood ElMaestro 2020-08-05 08:13
- Not understood PharmCat 2020-08-05 16:37
- Open issues ElMaestro 2020-08-06 21:11
- Open issues PharmCat 2020-08-07 00:02
- Open issues ElMaestro 2020-08-07 07:49
- Open issues PharmCat 2020-08-07 11:29
- Still can't make it work ElMaestro 2020-08-07 13:08
- Still can't make it work PharmCat 2020-08-07 15:42
- Still can't make it work ElMaestro 2020-08-07 16:44
- Still can't make it workPharmCat 2020-08-07 18:14
- Still can't make it work ElMaestro 2020-08-07 18:23
- And now it works ElMaestro 2020-08-07 21:31
- Still can't make it work PharmCat 2020-08-08 01:08
- Speed improvement ElMaestro 2020-08-08 12:25
- Speed improvement PharmCat 2020-08-08 17:27
- Speed improvement ElMaestro 2020-08-08 18:10
- Speed improvement PharmCat 2020-08-09 18:22
- Some tests... PharmCat 2020-08-10 11:48
- Speed improvement ElMaestro 2020-08-08 12:25
- Still can't make it work ElMaestro 2020-08-07 18:23
- Still can't make it workPharmCat 2020-08-07 18:14
- Still can't make it work ElMaestro 2020-08-07 16:44
- Still can't make it work PharmCat 2020-08-07 15:42
- Still can't make it work ElMaestro 2020-08-07 13:08
- Open issues PharmCat 2020-08-07 11:29
- Open issues ElMaestro 2020-08-07 07:49
- Open issues PharmCat 2020-08-07 00:02
- Open issues ElMaestro 2020-08-06 21:11
- Not understood PharmCat 2020-08-05 16:37
- Not understood ElMaestro 2020-08-05 08:13
- Not understood PharmCat 2020-08-05 01:41
- Not understood ElMaestro 2020-08-03 22:55
- And by the way.... PharmCat 2020-08-03 14:24
- And by the way.... ElMaestro 2020-07-24 12:52
- Semireplicated + REML in R Helmut 2020-07-10 19:03