Skip to content

Instantly share code, notes, and snippets.

@bbolker
Created July 16, 2020 15:48
Show Gist options
  • Save bbolker/501c1dc6b1ed8be5c46ef66da5a964aa to your computer and use it in GitHub Desktop.
Save bbolker/501c1dc6b1ed8be5c46ef66da5a964aa to your computer and use it in GitHub Desktop.
library(nlme)
packageVersion('nlme')
## the values are fixed for all models
rho <- 0.3
fixform <- distance ~ age + factor(Sex)
## remove the grouping and other classes to keep it simple
## full balanced data set
Obal<-as.data.frame(Orthodont)
## create unbalanced Orthodont data set with no singletons
Ounbal <-as.data.frame(Orthodont[c(3,4,5, 7,8, 11, 12, 15, 16, 19, 20:28,31,32,35,36,39,
40,41,43,44,47,48,51,52,53,55,56,59,60,63,64,67,68,71,72,75,76,
79,80,83:85,87,88,91,92,95,96,99:101,103:105,107,108),])
## sparser unbalanced data set with singletons
Ounbal2<-as.data.frame(Orthodont[c(1,5,9,13,17,18,21:23,30:31,44,45:47,56:58,65:71,74:76,77,81,85,89,90,
93:94, 97,101:108),])
Subjs<-c(paste0(c("M"),sprintf("%02d",1:16)), paste0(c("F"),sprintf("%02d",1:11)))
table(Obal$Subject)[Subjs]
table(Ounbal$Subject)[Subjs]
table(Ounbal2$Subject)[Subjs]
## test with the full balanced data set
## pre-calling Initialize doesn't seem to make a difference
csB <- corCompSymm(value = rho, form = ~ 1 | Subject)
csB<-Initialize(csB, data=Obal)
print(mod_bal_init<-gls(fixform, correlation=csB, data=Obal))
## compare with direct call
print(mod_bal_noinit <-gls(fixform, correlation=corCompSymm(value = rho, form = ~ 1 | Subject), data=Obal) )
##--------------
## compare Unbalanced:
csU <- corCompSymm(value = rho, form = ~ 1 | Subject)
csU<-Initialize(csU, data=Ounbal)
print(mod_ubal1_init<-gls(fixform, correlation=csU, data=Ounbal) )
## compare with direct call
print(mod_ubal1_noinit<-gls(fixform, correlation=corCompSymm(value = rho, form = ~ 1 | Subject), data=Ounbal) )
## compare with Unbalanced 2 with singletons
csU2 <- corCompSymm(value = rho, form = ~ 1 | Subject)
csU2<-Initialize(csU2, data=Ounbal2)
print(mod_ubal2_init<-gls(fixform, correlation=csU2, data=Ounbal2))
## compare with direct call
print(mod_ubal2_noinit<-gls(fixform, correlation=corCompSymm(value = rho, form = ~ 1 | Subject), data=Ounbal2))
mod_bal_lme <- lme(fixform, random=~1|Subject, data=Obal)
mod_ubal1_lme <- lme(fixform, random=~1|Subject, data=Ounbal)
mod_ubal2_lme <- lme(fixform, random=~1|Subject, data=Ounbal2)
mnames <- ls(pattern="^mod_")
mList <- mget(mnames)
library(broom.mixed)
library(tidyverse)
comb <- (purrr::map_dfr(mList,tidy,.id="model",conf.int=TRUE)
%>% filter(term=="age")
%>% select(model,estimate,std.error,conf.low,conf.high)
%>% separate(model,into=c("junk","data","method"))
%>% select(-junk)
%>% mutate_at("method",factor,levels=c("init","noinit","lme"))
)
library(ggplot2); theme_set(theme_bw())
library(ggstance)
library(colorspace)
scale_colour_discrete <- scale_colour_discrete_qualitative
ggplot(comb,aes(x=estimate,xmin=conf.low,xmax=conf.high,y=method,
colour=method)) +
geom_pointrange(position=position_dodgev(height=0.5))+
facet_wrap(~data,ncol=1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment