Merge pull request #49 from PharmCat/dev

repeated UN warn

Project.toml

@@ -3,7 +3,7 @@ uuid = "a1dec852-9fe5-11e9-361f-8d9fde67cfa2"
 keywords = ["lenearmodel", "mixedmodel"]
 desc = "Mixed-effects models with flexible covariance structure."
 authors = ["Vladimir Arnautov <[email protected]>"]
-version = "0.16.0"
+version = "0.16.1"
 
 [deps]
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"

docs/src/api.md

@@ -124,12 +124,6 @@ Metida.ToeplitzParameterized
 Metida.Unstructured
 ```
 
-### Metida.ScaledWeightedCov
-
-```@docs
-Metida.ScaledWeightedCov
-```
-
 ### Methods
 
 ### Metida.caic
@@ -388,6 +382,12 @@ Metida.SpatialGaussianD
 Metida.SpatialPowerD
 ```
 
+### Metida.ScaledWeightedCov
+
+```@docs
+Metida.ScaledWeightedCov
+```
+
 ### Metida.dof_contain
 
 ```@docs

src/lmm.jl

@@ -84,7 +84,9 @@ struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
             union!(tv, (wts,))
         end
         ct = Tables.columntable(data)
-        if !(tv ⊆ keys(ct)) error("Some column(s) not found!") end
+        if !(tv ⊆ keys(ct)) 
+            error("Column(s) ($(setdiff(tv, keys(ct))) not found in table!") 
+        end
         data, data_ = StatsModels.missing_omit(NamedTuple{tuple(tv...)}(ct))
         lmmlog = Vector{LMMLogMsg}(undef, 0)
         sch    = schema(model, data, contrasts)
@@ -151,10 +153,10 @@ struct LMM{T <: AbstractFloat, W <: Union{LMMWts, Nothing}} <: MetidaModel
 
         mres = ModelResult(false, nothing, fill(NaN, covstr.tl), NaN, fill(NaN, coefn), nothing, fill(NaN, coefn, coefn), fill(NaN, coefn), nothing, false)
 
-        LMM(model, f, ModelStructure(assign), covstr, lmmdata, LMMDataViews(lmmdata.xv, lmmdata.yv, covstr.vcovblock), nfixed, rankx, mres, findmax(length, covstr.vcovblock)[1], lmmwts, lmmlog)
+        return LMM(model, f, ModelStructure(assign), covstr, lmmdata, LMMDataViews(lmmdata.xv, lmmdata.yv, covstr.vcovblock), nfixed, rankx, mres, findmax(length, covstr.vcovblock)[1], lmmwts, lmmlog)
     end
     function LMM(f::LMMformula, data; kwargs...)
-        LMM(f.formula, data; random = f.random, repeated = f.repeated, kwargs...)
+        return LMM(f.formula, data; random = f.random, repeated = f.repeated, kwargs...)
     end
 end
 
@@ -165,7 +167,7 @@ end
 Length of theta vector.
 """
 function thetalength(lmm)
-    lmm.covstr.tl
+    return lmm.covstr.tl
 end
 
 """
@@ -174,7 +176,7 @@ end
 Coef number.
 """
 function coefn(lmm)
-    length(lmm.result.beta)
+    return length(lmm.result.beta)
 end
 
 """
@@ -183,28 +185,28 @@ end
 Return theta vector.
 """
 function theta(lmm::LMM)
-    copy(theta_(lmm))
+    return copy(theta_(lmm))
 end
 function theta_(lmm::LMM)
-    lmm.result.theta
+    return lmm.result.theta
 end
 """
     rankx(lmm::LMM)
 
 Return rank of `X` matrix.
 """
 function rankx(lmm::LMM)
-    Int(lmm.rankx)
+    return Int(lmm.rankx)
 end
 
 function nblocks(mm::MetidaModel)
     return length(mm.covstr.vcovblock)
 end
 function maxblocksize(mm::MetidaModel)
-    mm.maxvcbl
+    return mm.maxvcbl
 end
 function assign(lmm::LMM)
-    lmm.modstr.assign
+    return lmm.modstr.assign
 end
 
 ################################################################################
@@ -220,17 +222,17 @@ function lmmlog!(io, lmmlog::Vector{LMMLogMsg}, verbose, vmsg)
     end
 end
 function lmmlog!(lmmlog::Vector{LMMLogMsg}, verbose, vmsg)
-    lmmlog!(stdout, lmmlog, verbose, vmsg)
+    return lmmlog!(stdout, lmmlog, verbose, vmsg)
 end
 function lmmlog!(io, lmm::LMM, verbose, vmsg)
-    lmmlog!(io, lmm.log, verbose, vmsg)
+    return lmmlog!(io, lmm.log, verbose, vmsg)
 end
 #MetidaNLopt use this
 function lmmlog!(lmm::LMM, verbose, vmsg)
-    lmmlog!(stdout, lmm, verbose, vmsg)
+    return lmmlog!(stdout, lmm, verbose, vmsg)
 end
 function lmmlog!(lmm::LMM, vmsg)
-    lmmlog!(stdout, lmm, 1, vmsg)
+    return lmmlog!(stdout, lmm, 1, vmsg)
 end
 
 function msgnum(log::Vector{LMMLogMsg}, type::Symbol)
@@ -240,10 +242,10 @@ function msgnum(log::Vector{LMMLogMsg}, type::Symbol)
             n += 1
         end
     end
-    n
+    return n
 end
 function msgnum(log::Vector{LMMLogMsg})
-    length(log)
+    return length(log)
 end
 ################################################################################
 
@@ -321,7 +323,7 @@ function Base.show(io::IO, lmm::LMM)
 end
 
 function printresult(io, res::T) where T <: Optim.MultivariateOptimizationResults
-    Optim.converged(res) ? printstyled(io, "converged"; color = :green) : printstyled(io, "not converged"; color = :red)
+    return Optim.converged(res) ? printstyled(io, "converged"; color = :green) : printstyled(io, "not converged"; color = :red)
 end
 function printresult(io, res)
     if res[3] == :FTOL_REACHED || res[3] == :XTOL_REACHED || res[3] == :SUCCESS
@@ -349,7 +351,7 @@ end
 Return fitting log.
 """
 function getlog(lmm::LMM)
-    lmm.log
+    return lmm.log
 end
 
 ################################################################################
@@ -360,7 +362,7 @@ function Base.getproperty(x::LMM, s::Symbol)
     elseif s == :β
         return x.result.beta
     end
-    getfield(x, s)
+    return getfield(x, s)
 end
 
 #=

src/varstruct.jl

@@ -316,6 +316,19 @@ struct CovStructure{T, T2} <: AbstractCovarianceStructure
             else
                 dicts[rn+i]  = Dict(1 => collect(1:rown)) #changed to range
             end
+            # If UN structure used all repeated levels should be unique within one subject, otherwise results can be meaningless!
+            wflag = true
+            if isa(repeated[i].covtype.s, UN_)
+                for (k,v) in dicts[rn+i]
+                    sv = view(rz_[i], v, :)
+                    for j = 1:size(sv, 2)
+                        if sum(view(sv, :, j)) > 1 && wflag
+                            wflag = false
+                            @warn "If UN structure used for repeated effect all levels should be unique within one subject, otherwise results can be meaningless!"
+                        end
+                    end
+                end
+            end
 
             sn[rn + i]   = length(dicts[rn+i])
             q[rn + i]    = size(rz_[i], 2)

src/vartypes.jl

@@ -389,7 +389,8 @@ const SPGAUD = SpatialGaussianD()
 """
     Unstructured()
 
-Unstructured covariance structure with `t*(t+1)/2-t` paremeters where `t` - number of factor levels, `t*(t+1)/2-2t` of them is covariance (ρ) patemeters.
+Unstructured covariance structure with `t*(t+1)/2-t` paremeters where `t` - number of factor levels, `t*(t+1)/2-2t` of them is covariance (ρ) patemeters. 
+All levels for repeated effect should be unique within each subject.
 
 UN = Unstructured()
 """

test/test.jl

@@ -609,10 +609,8 @@ end
     @test Metida.m2logreml(lmm)  ≈ 713.5850978377632 atol=1E-8
 end
 @testset "  Model: BE RDS 1, FDA model                               "  begin
-    dfrds        = CSV.File(joinpath(path, "csv", "berds", "rds1.csv"), types = Dict(:PK => Float64, :subject => String, :period => String, :sequence => String, :treatment => String )) |> DataFrame
-    dropmissing!(dfrds)
-    dfrds.lnpk = log.(dfrds.PK)
-    lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrds;
+
+    lmm = Metida.LMM(@formula(lnpk~sequence+period+treatment), dfrdsfda;
     random = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.CSH),
     repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.DIAG),
     )
@@ -627,7 +625,7 @@ end
     @test est.cil[1] ≈ 0.06863 atol=1E-4
     @test est.ciu[1] ≈ 0.2223 atol=1E-4
 
-    lmm = Metida.LMM(@formula(lnpk~0+sequence+period+treatment), dfrds;
+    lmm = Metida.LMM(@formula(lnpk~0+sequence+period+treatment), dfrdsfda;
     random = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.CSH),
     repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.DIAG),
     )
@@ -651,6 +649,15 @@ end
     repeated = Metida.VarEffect(Metida.@covstr(Formulation|Subject), Metida.UN),
     )
     Metida.fit!(lmm)
+    @test Metida.m2logreml(lmm)  ≈ -3.895979534278979 atol=1E-8
+
+
+    lmm2 = Metida.LMM(@formula(log(Var)~Sequence+Period+Formulation), dfrds;
+    repeated = Metida.VarEffect(Metida.@covstr(Formulation|Subject), Metida.CSH),
+    )
+    Metida.fit!(lmm2)
+
+    @test Metida.m2logreml(lmm)  ≈ Metida.m2logreml(lmm2) atol=1E-8
 end
 
 
@@ -845,6 +852,10 @@ end
     )
     Metida.fit!(lmm)
     println(io, lmm.log)
+
+    # Warn for non-unique levels for repeated effect within subject
+    @test_warn "If UN structure used for repeated effect all levels should be unique within one subject, otherwise results can be meaningless!" Metida.LMM(@formula(log(lnpk)~sequence+period+treatment), dfrdsfda; repeated = Metida.VarEffect(Metida.@covstr(treatment|subject), Metida.UN))
+
 end
 ################################################################################
 #                                  Sweep test

test/testdata.jl

@@ -16,3 +16,8 @@ ftdf3        = CSV.File(path*"/csv/ftdf3.csv"; types =
 [String,  Float64, Float64, String, String, String, String, String, Float64]) |> DataFrame
 
 spatdf       = CSV.File(path*"/csv/spatialdata.csv"; types = [Int, Int, String, Float64, Float64, Float64, Float64, Float64]) |> DataFrame
+
+
+dfrdsfda        = CSV.File(joinpath(path, "csv", "berds", "rds1.csv"), types = Dict(:PK => Float64, :subject => String, :period => String, :sequence => String, :treatment => String )) |> DataFrame
+dropmissing!(dfrdsfda)
+dfrdsfda.lnpk = log.(dfrdsfda.PK)