Skip to content

Commit

Permalink
update, no random
Browse files Browse the repository at this point in the history
  • Loading branch information
PharmCat committed Jan 6, 2021
1 parent d80d56b commit d8663b7
Show file tree
Hide file tree
Showing 7 changed files with 73 additions and 26 deletions.
15 changes: 13 additions & 2 deletions src/gmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,8 @@ function gmat_switch!(G, θ, covstr, i)
gmat_csh!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype)
elseif covstr.random[i].covtype.s == :CS
gmat_cs!(G, θ[covstr.tr[i]], covstr.q[i], covstr.random[i].covtype)
elseif covstr.random[i].covtype.s == :ZERO
gmat_zero!(G, similar(θ, 0), covstr.q[i], covstr.random[i].covtype)
else
throw(ErrorException("Unknown covariance structure: $(covstr.random[i].covtype.s), n = $(i)"))
end
Expand Down Expand Up @@ -76,6 +78,10 @@ end


################################################################################
function gmat_zero!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T
mx .= zero(T)
nothing
end
function gmat_si!(mx, θ::Vector{T}, zn::Int, ::CovarianceType) where T
val = θ[1] ^ 2
for i = 1:size(mx, 1)
Expand All @@ -90,11 +96,16 @@ function gmat_diag!(mx, θ::Vector{T}, ::Int, ::CovarianceType) where T
nothing
end
function gmat_ar!(mx, θ::Vector{T}, zn::Int, ::CovarianceType) where T
mx .= θ[1] ^ 2
#mx .= θ[1] ^ 2
de = θ[1] ^ 2
for i = 1:zn
mx[i, i] = de
end
if zn > 1
for m = 1:zn - 1
for n = m + 1:zn
@inbounds mx[m, n] = mx[m, m] * θ[2] ^ (n - m)
ode = de * θ[2] ^ (n - m)
@inbounds mx[m, n] = ode
@inbounds mx[n, m] = mx[m, n]
end
end
Expand Down
9 changes: 6 additions & 3 deletions src/lmm.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ struct LMM{T} <: MetidaModel
mm = ModelMatrix(mf)
if random === nothing
random = VarEffect()
#random = RZero()
end
if repeated === nothing
repeated = VarEffect()
Expand Down Expand Up @@ -80,12 +81,12 @@ function Base.show(io::IO, lmm::LMM)
println(io, "Random $i: ")
println(io, " Model: $(lmm.covstr.random[i].model === nothing ? "nothing" : lmm.covstr.random[i].model)")
println(io, " Type: $(lmm.covstr.random[i].covtype.s) ($(lmm.covstr.t[i]))")
println(io, " Coefnames: $(coefnames(lmm.covstr.schema[i]))")
#println(io, " Coefnames: $(coefnames(lmm.covstr.schema[i]))")
end
println(io, "Repeated: ")
println(io, " Model: $(lmm.covstr.repeated.model === nothing ? "nothing" : lmm.covstr.repeated.model)")
println(io, " Type: $(lmm.covstr.repeated.covtype.s) ($(lmm.covstr.t[end]))")
println(io, " Coefnames: $(lmm.covstr.repeated.model === nothing ? "-" : coefnames(lmm.covstr.schema[end]))")
#println(io, " Coefnames: $(lmm.covstr.repeated.model === nothing ? "-" : coefnames(lmm.covstr.schema[end]))")
println(io, "")
if lmm.result.fit
print(io, "Status: ")
Expand All @@ -110,7 +111,9 @@ function Base.show(io::IO, lmm::LMM)
mx = hcat(Matrix{Any}(undef, lmm.covstr.tl, 1), lmm.covstr.rcnames, lmm.covstr.ct, round.(lmm.result.theta, sigdigits = 6))

for i = 1:length(lmm.covstr.random)
view(mx, lmm.covstr.tr[i], 1) .= "Random $i"
if lmm.covstr.random[i].covtype.s != :ZERO
view(mx, lmm.covstr.tr[i], 1) .= "Random $i"
end
end
view(mx, lmm.covstr.tr[end], 1) .= "Repeated"
for i = 1:lmm.covstr.tl
Expand Down
2 changes: 1 addition & 1 deletion src/rmat.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ function rmatp_diag!(mx, θ::Vector{T}, rz, ::CovarianceType) where T
end
nothing
end
function rmatp_ar!(mx, θ::Vector{T}, rz, ::CovarianceType) where T
function rmatp_ar!(mx, θ::Vector{T}, ::AbstractMatrix, ::CovarianceType) where T
rn = size(mx, 1)
de = θ[1] ^ 2
for m = 1:rn
Expand Down
32 changes: 24 additions & 8 deletions src/varstruct.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,10 @@ function HeterogeneousCompoundSymmetry()
end
const CSH = HeterogeneousCompoundSymmetry()

function RZero()
CovarianceType(:ZERO, x -> 0, x -> 0, x -> 0)
end

#TOE

#TOEH
Expand Down Expand Up @@ -114,8 +118,8 @@ struct VarEffect
function VarEffect(model; coding = nothing)
VarEffect(model, DIAG, coding)
end
function VarEffect(covtype::T; coding = nothing) where T <: AbstractCovarianceType
VarEffect(@covstr(1), covtype, coding)
function VarEffect(covtype::T; coding = nothing, subj = nothing) where T <: AbstractCovarianceType
VarEffect(@covstr(1), covtype, coding; subj = subj)
end
function VarEffect()
VarEffect(@covstr(1), SI, Dict{Symbol, AbstractContrasts}())
Expand Down Expand Up @@ -170,6 +174,13 @@ struct CovStructure{T} <: AbstractCovarianceStructure
#
# RANDOM EFFECTS
for i = 1:length(random)
if random[i].covtype.s == :ZERO
q[i] = 0
t[i] = 0
tr[i] = 0:0
zrndur[i] = 0:0
continue
end
if length(random[i].coding) == 0 && random[i].fulldummy
fill_coding_dict!(random[i].model, random[i].coding, data)
end
Expand Down Expand Up @@ -206,6 +217,9 @@ struct CovStructure{T} <: AbstractCovarianceStructure
rcnames = Vector{String}(undef, tl)
ctn = 1
for i = 1:length(random)
if random[i].covtype.s == :ZERO
continue
end
for i2 = 1:random[i].covtype.v(q[i])
ct[ctn] = :var
ctn +=1
Expand All @@ -229,12 +243,14 @@ struct CovStructure{T} <: AbstractCovarianceStructure
end
end
view(rcnames, tr[end]) .= rcoefnames(schema[end], t[end], Val{repeated.covtype.s}())
for i = 1:length(blocks)
sblock[i] = Vector{Vector{Vector{UInt32}}}(undef, alleffl)
for s = 1:alleffl
sblock[i][s] = Vector{Vector{UInt32}}(undef, 0)
for col in eachcol(view(subjz[s], blocks[i], :))
if any(col) push!(sblock[i][s], findall(x->x==true, col)) end
if random[1].covtype.s != :ZERO
for i = 1:length(blocks)
sblock[i] = Vector{Vector{Vector{UInt32}}}(undef, alleffl)
for s = 1:alleffl
sblock[i][s] = Vector{Vector{UInt32}}(undef, 0)
for col in eachcol(view(subjz[s], blocks[i], :))
if any(col) push!(sblock[i][s], findall(x->x==true, col)) end
end
end
end
end
Expand Down
26 changes: 23 additions & 3 deletions test/ar.jl
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,33 @@


#=
SPSS
REML 453.339544
=#
df = CSV.File(path*"/csv/RepeatedPulse.csv") |> DataFrame
df.Pulse = float.(df.Pulse)
categorical!(df, :Time);
categorical!(df, :Day);
@testset " RepeatedPulse.csv" begin
lmm = Metida.LMM(@formula(Pulse~1), df;
random = Metida.VarEffect(Metida.@covstr(Time), Metida.SI),
repeated = Metida.VarEffect(Metida.@covstr(1), Metida.AR),
repeated = Metida.VarEffect(Metida.@covstr(Day), Metida.AR),
subject = :Time
)
Metida.fit!(lmm)
@test lmm.result.reml ≈ 450.75260020816347 atol=1E-6
end

lmm = Metida.LMM(@formula(Pulse~1), df;
random = Metida.VarEffect(Metida.@covstr(1), Metida.AR),
subject = :Time
)
Metida.fit!(lmm)
@test lmm.result.reml ≈ 453.3395546350102 atol=1E-6
end

#=
SPSS
REML 5456.638120
=#
df = CSV.File(path*"/csv/ChickWeight.csv") |> DataFrame
df.weight = float.(df.weight)
df.Time = float.(df.Time)
Expand All @@ -30,4 +43,11 @@ categorical!(df, :Diet);
)
Metida.fit!(lmm)
@test lmm.result.reml ≈ 5453.137479692881 atol=1E-6

lmm = Metida.LMM(@formula(weight~1 + Diet & Time), df;
random = Metida.VarEffect(Metida.@covstr(Diet), Metida.ARH),
subject = :Chick
)
Metida.fit!(lmm)
@test lmm.result.reml ≈ 5480.751465914058 atol=1E-6
end
12 changes: 5 additions & 7 deletions test/lme4.jl
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,16 @@


df = CSV.File(path*"/csv/lme4/sleepstudy.csv") |> DataFrame

categorical!(df, :Subject);
categorical!(df, :Days);


#=
SPSS
REML 1729.492560
=#
@testset " sleepstudy.csv" begin
lmm = Metida.LMM(@formula(Reaction~Days), df;
random = Metida.VarEffect(Metida.@covstr(1), Metida.SI),
random = Metida.VarEffect(Metida.SI),
subject = :Subject
)
Metida.fit!(lmm)
Expand All @@ -26,10 +24,10 @@ categorical!(df, :sample);

@testset " Penicillin.csv" begin
lmm = Metida.LMM(@formula(diameter~1), df;
random = [Metida.VarEffect(Metida.@covstr(1), Metida.SI, subj = :plate), Metida.VarEffect(Metida.@covstr(1), Metida.SI, subj = :sample)]
random = [Metida.VarEffect(Metida.SI, subj = :plate), Metida.VarEffect(Metida.SI, subj = :sample)]
)
Metida.fit!(lmm)
#@test lmm.result.reml ≈ 330.86058899109184 atol=1E-6`
@test lmm.result.reml ≈ 330.86058899109184 atol=1E-6
end


Expand All @@ -39,8 +37,8 @@ categorical!(df, :cask);
categorical!(df, :sample);
@testset " Pastes.csv" begin
lmm = Metida.LMM(@formula(strength~1), df;
random = [Metida.VarEffect(Metida.@covstr(1), Metida.SI, subj = :batch), Metida.VarEffect(Metida.@covstr(1), Metida.SI, subj = [:batch, :cask])]
random = [Metida.VarEffect(Metida.SI, subj = :batch), Metida.VarEffect(Metida.SI, subj = [:batch, :cask])]
)
Metida.fit!(lmm)
#@test lmm.result.reml ≈ 246.99074585348623 atol=1E-6`
@test lmm.result.reml ≈ 246.99074585348623 atol=1E-6
end
3 changes: 1 addition & 2 deletions test/test.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,7 +78,6 @@ end
)
Metida.fit!(lmm)
@test Metida.m2logreml(lmm) ≈ 10.065239006121315 atol=1E-6

end
@testset " Model: SI/subject + DIAG " begin
lmm = Metida.LMM(@formula(var~sequence+period+formulation), df0;
Expand Down Expand Up @@ -126,7 +125,7 @@ end
random = Metida.VarEffect(Metida.@covstr(formulation * period), Metida.DIAG), subject = :subject
)
Metida.fit!(lmm)
#@test lmm.result.reml ≈ 13.555817544390917 atol=1E-4
#@test lmm.result.reml ≈ 13.555817544390917 atol=1E-4
@test true
end

Expand Down

0 comments on commit d8663b7

Please sign in to comment.