Skip to content

Commit

Permalink
Merge pull request #10 from yufongpeng/release-0.2.4
Browse files Browse the repository at this point in the history
Release 0.2.4
  • Loading branch information
yufongpeng committed Feb 28, 2024
2 parents 9293537 + d0d43a5 commit fd7e625
Show file tree
Hide file tree
Showing 5 changed files with 98 additions and 13 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/CI.yml
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ jobs:
fail-fast: false
matrix:
version:
- '1.6'
- '1.8'
- '1'
- 'nightly'
os:
Expand Down
6 changes: 3 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand All @@ -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"
14 changes: 7 additions & 7 deletions src/anova.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

#=
Expand Down Expand Up @@ -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} =
Expand Down Expand Up @@ -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}} =
Expand Down
4 changes: 2 additions & 2 deletions src/fit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)]]
Expand Down Expand Up @@ -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;
Expand Down
85 changes: 85 additions & 0 deletions test/calcDenDF.R
Original file line number Diff line number Diff line change
@@ -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
}

2 comments on commit fd7e625

@yufongpeng
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/101930

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.2.4 -m "<description of version>" fd7e62552524e5e369fe029464bf7102bcdd45c0
git push origin v0.2.4

Please sign in to comment.