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

posted by PharmCat  – Russia, 2020-08-07 16:14 (169 d 23:56 ago) – Posting: # 21835
Views: 6,820

(edited by PharmCat on 2020-08-07 23:12)

» 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:

Activity
 Admin contact
21,310 posts in 4,445 threads, 1,489 registered users;
online 3 (0 registered, 3 guests [including 1 identified bots]).
Forum time: Sunday 16:11 UTC (Europe/Vienna)

Every man gets a narrower and narrower field of knowledge
in which he must be an expert in order to compete with other people.
The specialist knows more and more about less and less
and finally knows everything about nothing.    Konrad Lorenz

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5