diff --git a/Project.toml b/Project.toml index 976f2a7..98ea8f3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "JellyMe4" uuid = "19ac8677-9a15-4623-9afd-84acc6165ce7" authors = ["Phillip Alday "] -version = "0.2.7" +version = "0.2.8" [deps] CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" @@ -17,16 +17,16 @@ StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [compat] -CategoricalArrays = "0.7, 0.8, 0.9, 0.10" -DataFrames = "0.20, 0.21, 0.22, 1.0" -Distributions = "0.20, 0.22, 0.23, 0.24, 0.25" +CategoricalArrays = "0.10" +DataFrames = "1" +Distributions = "0.25" GLM = "1.3" -MixedModels = "3.0, 4.0" +MixedModels = "4" RCall = "0.13" RegressionFormulae = "0.1" StatsBase = "0.31, 0.32, 0.33" -StatsModels = "0.6.8" -Tables = "0.2, 1.0" +StatsModels = "0.6, 0.7" +Tables = "1" julia = "1.6" [extras] diff --git a/deps/build.jl b/deps/build.jl index 270776e..5be2147 100644 --- a/deps/build.jl +++ b/deps/build.jl @@ -1 +1 @@ -# should we try to install lme4 or even afex automatically at build time? \ No newline at end of file +# should we try to install lme4 or even afex automatically at build time? diff --git a/src/glmerMod.jl b/src/glmerMod.jl index 7554e17..b6e9b19 100644 --- a/src/glmerMod.jl +++ b/src/glmerMod.jl @@ -119,7 +119,7 @@ function rcopy(::Type{GeneralizedLinearMixedModel}, s::Ptr{S4Sxp}) if length(θ) != length(m.θ) @error """You're probably using || in R with a categorical variable, - whose translation is currently unsupported with MixedModels 3.0.""" + whose translation is currently unsupported with recent MixedModels releases""" throw(ArgumentError("Parameter vectors in R and Julia are different sizes.")) end diff --git a/src/lmerMod.jl b/src/lmerMod.jl index 8e97fd8..2846480 100644 --- a/src/lmerMod.jl +++ b/src/lmerMod.jl @@ -57,7 +57,7 @@ function rcopy(::Type{LinearMixedModel}, s::Ptr{S4Sxp}) if length(θ) != length(m.θ) @error """You're probably using || in R with a categorical variable, - whose translation is currently unsupported with MixedModels 3.0.""" + whose translation is currently unsupported with recent MixedModels releases""" throw(ArgumentError("Parameter vectors in R and Julia are different sizes.")) end diff --git a/src/utilities.jl b/src/utilities.jl index 6155389..756f906 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -41,16 +41,30 @@ function get_r_contrasts(rdf) data = rcopy(rdf) # get categorical columns cnames = [c for c in propertynames(data) if typeof(data[!, c]) <: CategoricalArray] - return Dict(c => HypothesisCoding(pinv(rcopy(R"contrasts($(rdf[c]))")); - labels=rcopyarray(R"colnames(contrasts($(rdf[c])))")) - for c in cnames) + pairs = [] + for c in cnames + @debug "" c + levels = rcopyarray(R"rownames(contrasts($(rdf[c])))") + @debug "" levels + if levels == [nothing] + # for some contrasts, the rownames aren't defined + # we don't prefer this method because it's possible + # that it's not quite right when people are using really + # fancy contrasts and it can fail + levels = rcopyarray(R"levels($(rdf[c]))") + end + labels = rcopyarray(R"colnames(contrasts($(rdf[c])))") + @debug "" labels + push!(pairs, + c => HypothesisCoding(pinv(rcopy(R"contrasts($(rdf[c]))")); levels, labels)) + end + return Dict(pairs) end """ set the contrasts on an R dataframe this of course has a side effect: it changes the R dataframe """ - function set_r_contrasts!(rdfname, formula) fixefform = first(formula.rhs) diff --git a/test/glmerMod.jl b/test/glmerMod.jl index 72a5df8..d450b53 100644 --- a/test/glmerMod.jl +++ b/test/glmerMod.jl @@ -132,7 +132,7 @@ logistic(x) = 1 / (1 + exp(-x)) # we max out at 1 scalar RE for the nAGQ tests jlmm = GLMM(@formula(r2 ~ 1 + anger + gender + btype + situ + (1 | subj)), dat, Bernoulli()) - fit!(jlmm; fast=false, nAGQ=9) + fit!(jlmm; fast=false, nAGQ=9, progress=false) jm = (jlmm, dat) @rput jm # @test_warn Regex(".*categorical.*") @rput jm; @@ -147,7 +147,7 @@ logistic(x) = 1 / (1 + exp(-x)) # we max out at 1 scalar RE for the nAGQ tests jlmm = GLMM(@formula(r2 ~ 1 + anger + gender + btype + situ + (1 | subj)), dat, Bernoulli()) - fit!(jlmm) + fit!(jlmm; progress=false) jm = Tuple([jlmm, dat]) @rput jm @@ -159,7 +159,7 @@ logistic(x) = 1 / (1 + exp(-x)) @testset "columntable" begin jlmm = GLMM(@formula(r2 ~ 1 + anger + gender + btype + situ + (1 | subj)), dat, Bernoulli()) - fit!(jlmm) + fit!(jlmm; progress=false) jm = Tuple([jlmm, columntable(dat)]) @rput jm end @@ -170,7 +170,7 @@ logistic(x) = 1 / (1 + exp(-x)) dat = dataset(:cbpp) dat.rate = dat.incid ./ dat.hsz jlmm = fit(MixedModel, @formula(rate ~ 1 + period + (1 | herd)), - dat, Binomial(); wts=float(dat.hsz), fast=true) + dat, Binomial(); wts=float(dat.hsz), fast=true, progress=false) jm = (jlmm, dat) @rput jm @@ -192,7 +192,7 @@ logistic(x) = 1 / (1 + exp(-x)) dat, Poisson()) jm = (jlmm, dat) - fit!(jlmm; fast=true) # problems with this one in fast=false + fit!(jlmm; fast=true, progress=false) # problems with this one in fast=false @rput jm # @test_warn Regex(".*categorical.*") @rput jm; @test rcopy(R"""jm@devcomp$dims["nAGQ"]""") == 0 @@ -228,7 +228,7 @@ logistic(x) = 1 / (1 + exp(-x)) jlmm = fit(MixedModel, @formula(r2 ~ 1 + anger + gender + btype + situ + (1 | subj) + (1 | item)), - dat, Bernoulli(); + dat, Bernoulli(); progress=false, contrasts=Dict(:gender => EffectsCoding(), :btypes => EffectsCoding())) jm = (jlmm, dat) diff --git a/test/lmerMod.jl b/test/lmerMod.jl index 8a1d544..0924417 100644 --- a/test/lmerMod.jl +++ b/test/lmerMod.jl @@ -1,9 +1,10 @@ -using RCall, MixedModels, Test +using RCall, JellyMe4, MixedModels, Test using StatsBase: zscore using StatsModels: SeqDiffCoding using Tables: columntable -import JellyMe4: _set_lmer, _set_afex_installed +using JellyMe4: _set_lmer, _set_afex_installed +using MixedModels: dataset const LMM = LinearMixedModel const GLMM = GeneralizedLinearMixedModel @@ -25,7 +26,7 @@ const GLMM = GeneralizedLinearMixedModel # this is available in MixedModels.dataset(:sleepstudy) but with different # capitalization than in R - sleepstudy = rcopy(R"sleepstudy") + sleepstudy = rcopy(R"lme4::sleepstudy") kb07 = dataset(:kb07) @@ -33,14 +34,14 @@ const GLMM = GeneralizedLinearMixedModel ### from R ### jlmm = fit!(LMM(@formula(Reaction ~ 1 + Days + (1 + Days | Subject)), sleepstudy); - REML=false) - rlmm = rcopy(R"m <- lme4::lmer(Reaction ~ 1 + Days + (1 + Days|Subject),sleepstudy,REML=FALSE)") + REML=false, progress=false) + rlmm = rcopy(R"m <- lme4::lmer(Reaction ~ 1 + Days + (1 + Days|Subject),lme4::sleepstudy,REML=FALSE)") @test jlmm.θ ≈ rlmm.θ atol = 0.001 @test objective(jlmm) ≈ objective(rlmm) atol = 0.001 @test fixef(jlmm) ≈ fixef(rlmm) atol = 0.001 - jlmm = refit!(jlmm; REML=true) + jlmm = refit!(jlmm; REML=true, progress=false) rlmm = rcopy(R"update(m, REML=TRUE)") @test jlmm.θ ≈ rlmm.θ atol = 0.001 @@ -49,14 +50,14 @@ const GLMM = GeneralizedLinearMixedModel @testset "merModLmerTest" begin jlmm = fit!(LMM(@formula(Reaction ~ 1 + Days + (1 + Days | Subject)), - sleepstudy); REML=false) + sleepstudy); REML=false, progress=false) rlmm = rcopy(R"m <- lmerTest::lmer(Reaction ~ 1 + Days + (1 + Days|Subject),sleepstudy,REML=FALSE)") @test jlmm.θ ≈ rlmm.θ atol = 0.001 @test objective(jlmm) ≈ objective(rlmm) atol = 0.001 @test fixef(jlmm) ≈ fixef(rlmm) atol = 0.001 - jlmm = refit!(jlmm; REML=true) + jlmm = refit!(jlmm; REML=true, progress=false) rlmm = rcopy(R"update(m, REML=TRUE)") @test jlmm.θ ≈ rlmm.θ atol = 0.001 @@ -70,7 +71,7 @@ const GLMM = GeneralizedLinearMixedModel # the scalar-vector distinction for θ is missing in R @testset "scalar RE" begin jlmm = fit!(LMM(@formula(Reaction ~ 1 + Days + (1 | Subject)), sleepstudy); - REML=false) + REML=false, progress=false) rlmm = rcopy(R"m <- lme4::lmer(Reaction ~ 1 + Days + (1|Subject),sleepstudy,REML=FALSE)") @test jlmm.θ ≈ rlmm.θ atol = 0.001 @@ -91,7 +92,7 @@ const GLMM = GeneralizedLinearMixedModel @testset "contrasts" begin reval(""" - data(cake) + cake <- lme4::cake cake\$rr <- with(cake, replicate:recipe) """) rlmm = rcopy(R"fm1 <- lme4::lmer(angle ~ recipe * temperature + (1|rr), cake, REML= FALSE)") @@ -103,7 +104,7 @@ const GLMM = GeneralizedLinearMixedModel @testset "caret" begin reval(""" - data(cake) + cake <- lme4::cake cake\$rr <- with(cake, replicate:recipe) """) rlmm = rcopy(R"fm1 <- lme4::lmer(angle ~ (recipe + temperature)^2 + (1|rr), cake, REML= FALSE)") @@ -130,7 +131,7 @@ const GLMM = GeneralizedLinearMixedModel rlmm = rcopy(R"m <- lme4::lmer(Reaction ~ 1 + Days + (1 + Days||Subject),sleepstudy,REML=FALSE)") jlmm = fit!(LMM(@formula(Reaction ~ 1 + Days + zerocorr(1 + Days | Subject)), - sleepstudy); REML=false) + sleepstudy); REML=false, progress=false) # as a cheat for comparing the covariance matrices, we use PCA @test only(rlmm.rePCA) ≈ only(jlmm.rePCA) atol = 0.05 end @@ -151,7 +152,7 @@ const GLMM = GeneralizedLinearMixedModel @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) @test rcopy(R"REMLcrit(jm)") ≈ objective(jlmm) - refit!(jlmm; REML=false) + refit!(jlmm; REML=false, progress=false) jm = (jlmm, sleepstudy) @rput jm @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) @@ -168,14 +169,14 @@ const GLMM = GeneralizedLinearMixedModel R"m <- lme4::lmer(log10(Reaction) ~ 1 + log(Days2) + (1 + log(Days2)|Subject),sleepstudy,REML=FALSE)" jlmm = fit!(LMM(@formula(log10(Reaction) ~ 1 + log(Days2) + (1 + log(Days2) | Subject)), - sleepstudy); REML=false) + sleepstudy); REML=false, progress=false) jm = (jlmm, sleepstudy) @rput jm @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) @test rcopy(R"deviance(jm)") ≈ objective(jlmm) jlmm = fit!(LMM(@formula(Reaction ~ 1 + round(Days) + (1 | Subject)), - sleepstudy); REML=false) + sleepstudy); REML=false, progress=false) jm = (jlmm, sleepstudy) @test_throws ArgumentError (@rput jm) @@ -184,7 +185,7 @@ const GLMM = GeneralizedLinearMixedModel @testset "sorting by n BLUPs vs. n groups" begin jlmm = fit(LMM, @formula(rt_trunc ~ 1 + (1 | subj) + (1 + load + prec + spkr | item)), - kb07) + kb07; progress=false) jm = (jlmm, kb07) @rput jm @test rcopy(R"fitted(jm)") ≈ fitted(jlmm) @@ -202,7 +203,8 @@ const GLMM = GeneralizedLinearMixedModel cake """)) jlmm = fit(MixedModel, @formula(angle ~ recipe * temperature + (1 | rr)), - cake; REML=false, contrasts=Dict(:temperature => SeqDiffCoding())) + cake; REML=false, progress=false, + contrasts=Dict(:temperature => SeqDiffCoding())) jm = (jlmm, cake) @rput jm @test fixef(jlmm) ≈ rcopy(R"fixef(jm)") @@ -212,7 +214,7 @@ const GLMM = GeneralizedLinearMixedModel machines = rcopy(R"as.data.frame(nlme::Machines)") jlmm = fit(MixedModel, @formula(score ~ 1 + Machine + (1 + fulldummy(Machine) | Worker)), - machines) + machines; progress=false) rlmm = (jlmm, machines) @rput rlmm rlmmrepca = rcopy(R"summary(rePCA(rlmm))$Worker$importance[3,]") @@ -227,7 +229,7 @@ const GLMM = GeneralizedLinearMixedModel _set_afex_installed(false) jlmm = fit!(LMM(@formula(Reaction ~ 1 + Days + zerocorr(1 + Days | Subject)), - sleepstudy); REML=false) + sleepstudy); REML=false, progress=false) rlmm = (jlmm, sleepstudy) @rput rlmm @test rcopy(R"""!is(rlmm,"merModLmerTest")""") @@ -241,7 +243,7 @@ const GLMM = GeneralizedLinearMixedModel jlmm = fit(MixedModel, @formula(score ~ 1 + Machine + zerocorr(0 + Machine | Worker)), - machines) + machines; progress=false) @testset "afex pre-enabled" begin _set_lmer("afex::lmer_alt") @@ -279,7 +281,7 @@ const GLMM = GeneralizedLinearMixedModel slp[!, :obs] .= 1:nrow(slp) fm1 = fit(MixedModel, @formula(reaction ~ 1 + days + (1 | obs)), - slp; + slp; progress=false, contrasts=Dict(:obs => Grouping())) rfm1 = (fm1, slp) warning = r"""Warning: number of levels of each grouping factor must be < number of observations \(problems: obs\) diff --git a/test/merMod.jl b/test/merMod.jl index c6acf45..75fa5ea 100644 --- a/test/merMod.jl +++ b/test/merMod.jl @@ -6,9 +6,9 @@ const GLMM = GeneralizedLinearMixedModel @testset "merMod" begin # this is available in MixedModels.dataset(:sleepstudy) but with different # capitalization than in R - sleepstudy = rcopy(R"sleepstudy") + sleepstudy = rcopy(R"lme4::sleepstudy") jlmm = fit!(LMM(@formula(Reaction ~ 1 + round(Days) + (1 | Subject)), sleepstudy); - REML=false) + REML=false, progress=false) @testset "bare model" begin @test_throws ArgumentError (@rput jlmm) end