Release 0.2.4

.github/workflows/CI.yml

@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - '1.8'
           - '1'
           - 'nightly'
         os:

Project.toml

@@ -1,7 +1,7 @@
 name = "AnovaMixedModels"
 uuid = "6832ca0c-c2dd-4d0c-851a-521ca5db42f5"
 authors = ["Yu-Fong Peng <sciphypar@gmail.com>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 AnovaBase = "946dddda-6a23-4b48-8e70-8e60d9b8d680"
@@ -16,11 +16,11 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d"
 
 [compat]
-AnovaBase = "0.7"
+AnovaBase = "0.8"
 Distributions = "0.23, 0.24, 0.25"
 GLM = "1.9"
 MixedModels = "4"
 Reexport = "0.2, 1"
 StatsBase = "0.33, 0.34"
 StatsModels = "0.7"
-julia = "1.8, 1.9"
+julia = "1.8, 1.9, 1.10"

src/anova.jl

@@ -80,10 +80,10 @@ function anova(::Type{FTest},
         end
     elseif aovm.type == 2
         fstat = ntuple(last(fullasgn) - offset) do fix
-            select1 = sort!(collect(select_super_interaction(fullpred, fix + offset)))
-            select2 = setdiff(select1, fix + offset)
-            select1 = findall(in(select1), fullasgn)
-            select2 = findall(in(select2), fullasgn)
+            s1 = sort!(collect(select_super_interaction(fullpred, fix + offset)))
+            s2 = setdiff(s1, fix + offset)
+            select1 = findall(in(s1), fullasgn)
+            select2 = findall(in(s2), fullasgn)
             (β[select1]' * (varβ[select1, select1] \ β[select1]) - β[select2]' * (varβ[select2, select2] \ β[select2])) / df[fix] * adjust
         end
     else 
@@ -97,7 +97,7 @@ function anova(::Type{FTest},
     pvalue = ntuple(lastindex(fstat)) do id
         ccdf(FDist(df[id], dfr[id]), abs(fstat[id]))
     end
-    AnovaResult{FTest}(aovm, df, ntuple(x->NaN, length(fstat)), fstat, pvalue, (dof_residual = dfr,))
+    AnovaResult(aovm, FTest, df, ntuple(x->NaN, length(fstat)), fstat, pvalue, (dof_residual = dfr,))
 end
 
 #=
@@ -148,7 +148,7 @@ function anova(::Type{LikelihoodRatioTest},
     ord = sortperm(collect(df))
     df = df[ord]
     models = models[ord]
-    lrt_nested(NestedModels{M}(models), df, deviance.(models), 1)
+    lrt_nested(NestedModels(models), df, deviance.(models), 1)
 end
 
 anova(::Type{LikelihoodRatioTest}, aovm::NestedModels{M}) where {M <: MixedModel} = 
@@ -185,7 +185,7 @@ function anova(
     # isnested is not part of _iscomparable:  
     # isnested = true 
     dev = (_criterion(m0), deviance.(models[2:end])...)
-    lrt_nested(MixedAovModels{Union{M, T}}(models), df, dev, 1)
+    lrt_nested(MixedAovModels{Union{M, T}, length(models)}(models), df, dev, 1)
 end
 
 anova(::Type{LikelihoodRatioTest}, aovm::MixedAovModels{M}) where {M <:  Union{GLM_MODEL, MixedModel}} = 

src/fit.jl

@@ -64,7 +64,7 @@ function dof_residual_pred(model::LinearMixedModel)
 
     # Assign df to each factors based on the level
     assign = asgn(fixs)
-    has_intercept(fixs) && popfirst!(assign)
+    hasintercept(fixs) && popfirst!(assign)
     offset = first(assign) - 1
     dfr = ntuple(last(assign) - offset) do i
         dfperlv[level[findfirst(==(i + offset), assign)]]
@@ -125,7 +125,7 @@ function nestedmodels(model::M; null::Bool = true, kwargs...) where {M <: Linear
         fit!(lmm; REML, progress = true)
         lmm
     end
-    NestedModels{M}(models..., model)
+    NestedModels(models..., model)
 end
 
 nestedmodels(::Type{<: LinearMixedModel}, f::FormulaTerm, tbl; 

test/calcDenDF.R

@@ -0,0 +1,85 @@
+#' implementation of "between-within" df-guessing algorithm
+#' based on Pinheiro and Bates 2000 description, but tries
+#' to do better with random-slopes models
+
+#' @examples
+#' calcDenDF(~age,"Subject",nlme::Orthodont)
+#' calcDenDF(~age,data=nlme::Orthodont,random=~1|Subject)
+#' calcDenDF(~age,data=nlme::Orthodont,random=~age|Subject)
+#' ## off by 1
+#' calcDenDF(~poly(age,2),data=nlme::Orthodont,
+#'            random=~poly(age,2)|Subject)
+#' library(nlme)
+#' lmeDF <- function(formula=distance~age,random=~1|Subject) {
+#'     mod <- lme(formula,random,data=Orthodont)
+#'     aa <- anova(mod)
+#'     return(setNames(aa[,"denDF"],rownames(aa)))
+#' }
+#' lmeDF()
+#' lmeDF(random=~age|Subject) ## wrong!
+#' @param fixed fixed-effect formula
+#' @param grps (character) vector of grouping vectors
+#' @param random random-effect formula (can be specified instead of grps)
+#' @param n0 ?? number of parameters fitted at the population level
+calcDenDF <- function(fixed,grps,data,random=NULL,n0=1) {
+  ## n0 is a hack (still not matching some values properly)
+  ## try to construct random effects hierarchy
+  ## see whether effects vary within groups
+  testlevel1 <- function(grp,dmat) {
+    vv <- apply(dmat,2,
+                function(x)
+                  tapply(x,list(grp),function(x) length(unique(x))))
+    apply(vv>1,2,any)
+  }
+  ##
+  vlist <- list()
+  if (!is.null(random)) {
+    grps <- character(0)
+    if (!is.list(random)) random <- list(random)
+    lapply(random,
+           function(x) {
+             v <- deparse(x[[2]][[2]]) ## left side of bar
+             g <- deparse(x[[2]][[3]]) ## grouping factor (NB may be wrong for complex groupings)
+             v2 <- setdiff(colnames(model.matrix(reformulate(v),data=data)),"(Intercept)")
+             vlist[[g]] <<- c(vlist[[g]],v2)
+             grps <<- union(grps,g)  ## NB may not be correctly ordered for multiple groups
+           })
+    ## hack to combine
+    vnames <- unlist(Map(rep,as.list(names(vlist)),sapply(vlist,length)))
+    vlist <- setNames(unlist(vlist),vnames)
+  }
+  ## get grouping factors even if they are expressed as interactions
+  fdata <- lapply(data,factor)
+  gfacs <- data.frame(setNames(lapply(grps,
+                                      function(x) eval(parse(text=x),list2env(fdata))),grps),
+                      check.names=FALSE)
+
+  X <- model.matrix(fixed,data=data)
+  np <- ncol(X)
+  has.int <- "(Intercept)" %in% colnames(X)  ## FIXME: not tested
+  ## apply level-variation function to each grouping factor
+  lvary <- t(sapply(gfacs,testlevel1,dmat=X))
+  ## add intercept and observation levels
+  lvary <- rbind(pop=c(FALSE,rep(TRUE,np-1)),
+                 lvary,
+                 obs=rep(FALSE,np))
+  ##
+  for (i in seq_along(vlist)) {
+    w <- which(names(vlist[i])==rownames(lvary))
+    lvary[w,vlist[i]] <- FALSE
+  }
+  ## find index of first FALSE value 
+  lev <- apply(lvary,2,function(x) which(!x)[1])
+  p <- table(factor(lev,levels=seq(nrow(lvary))))  ## number of parameters estimated at each level
+  N <- nrow(X)
+  n <- c(pop=n0,sapply(gfacs,function(x) length(unique(x))),obs=N)  ## number of obs at each level
+  ndiff <- n[-1] - n[-length(n)]
+  dfs <- c(NA,setNames(ndiff-p[-1],names(p[-1])))
+  r <- setNames(dfs[as.character(lev)],names(lev))
+  if (has.int) {
+    r["(Intercept)"] <- dfs[length(dfs)]
+  } else {
+    r <- r[-1]
+  }
+  r
+}