[prev in list] [next in list] [prev in thread] [next in thread]
List: r-sig-mixed-models
Subject: Re: [R-sig-ME] Struggling (probably problem in wetware not software)
From: Phillip Alday <phillip.alday () mpi ! nl>
Date: 2020-05-03 23:11:27
Message-ID: 48545bce-b4b0-392c-b1af-7d6ac32ae669 () mpi ! nl
[Download RAW message or body]
It might make sense to have an item-level random effect, but you'll
probably need to exclude any covariates which are perfectly correlated
with a single item (e.g. risk item for item #4), because it will be very
hard to separate the covariate effect from the item effect.
Phillip
On 2/5/20 2:04 pm, Juho Kristian Ruohonen wrote:
> (reposting my response since the first version failed to include the
> mailing list among the recipients)
>
> Here's a non-statistician's take:
>
> You don't need a multilevel model for this. Given that every respondent
> contributes only one data point, no one's idiosyncracies are
> overrepresented, and we can expect them to cancel each other out. Your
> question can be modeled using standard logistic regression.
>
> Since you suspect *posItem* may influence the chance of omission, you
> should include it as an additional fixed effect besides *negItem*.
> Something like:
>
> mymodel <- glm(missing ~ posItem + riskItem, family = binomial, data =
> mydata)
>
> It'll be interesting to hear whether real statisticians disagree.
>
> Best,
>
> J
>
> pe 1. toukok. 2020 klo 21.38 Chris Evans (chrishold@psyctc.org) kirjoitti:
>
> > I am sorry if this is embarrassingly stupid/easy but I am failing to
> > understand something I thought would be fairly easy.
> >
> > The question of interest to me and others is whether people omit some
> > items of a questionnaire more than others.
> > To give a simple example I have several datasets for a ten item
> > self-report questionnaire where one item is about risk
> > and where three items are cued positively and the other seven are cued
> > negatively.
> >
> > Here is an example of some such data:
> >
> > ID Item Score missing missRand missReg riskItem posItem
> > <chr> <chr> <dbl> <dbl> <int> <int> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0 0 0 0 0
> > 2 ALCA0001 YPCORE02 0 0 0 0 0 0
> > 3 ALCA0001 YPCORE03 0 0 0 0 0 1
> > 4 ALCA0001 YPCORE04 0 0 0 0 1 0
> > 5 ALCA0001 YPCORE05 0 0 0 0 0 1
> > 6 ALCA0001 YPCORE06 1 0 0 0 0 0
> > 7 ALCA0001 YPCORE07 0 0 0 0 0 0
> > 8 ALCA0001 YPCORE08 0 0 0 0 0 0
> > 9 ALCA0001 YPCORE09 0 0 0 0 0 0
> > 10 ALCA0001 YPCORE10 0 0 0 0 0 1
> > 11 ALCA0004 YPCORE01 1 0 0 0 0 0
> > 12 ALCA0004 YPCORE02 0 0 0 0 0 0
> > 13 ALCA0004 YPCORE03 0 0 0 0 0 1
> > 14 ALCA0004 YPCORE04 0 0 0 0 1 0
> > 15 ALCA0004 YPCORE05 0 0 0 0 0 1
> >
> > missing is the observed omission of the item (1 = omitted)
> > missRand is a binomial random with probability = .5 I created
> > missReg is a biomial random with probability = item index number / 11,
> > i.e. strong relationship with item
> > riskItem identifies the risk items (item 4)
> > posItem identifies the positively cued items (3, 5 and 10)
> >
> > For this particular dataset there are 237 participants. (Like Groucho
> > Marks and principles, I have
> > others if you don't like this one!)
> >
> > The catch is that each participant has only completed the measure once.
> > It is a plausible model that
> > participants will vary in their propensity to omit items, it is also
> > highly plausible that participants
> > might omit the risk item more than the other item. It's a question of
> > interest (though I think not of
> > huge interest nor would I bet on it (!) that the positively cued items
> > might be omitted more than the
> > negatively cued ones or v.v.
> >
> > I thought I could test for these effects using glmer(). That works for
> >
> > > fitRisk <- glmer(missing ~ riskItem + (1 | ID), family = binomial, data
> > = longDat)
> > > summary(fitRisk)
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: missing ~ riskItem + (1 | ID)
> > Data: longDat
> >
> > AIC BIC logLik deviance df.resid
> > 51.7 69.0 -22.9 45.7 2367
> >
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -0.52020 -0.00113 -0.00113 -0.00113 2.93829
> >
> > Random effects:
> > Groups Name Variance Std.Dev.
> > ID (Intercept) 205.3 14.33
> > Number of obs: 2370, groups: ID, 237
> >
> > Fixed effects:
> > Estimate Std. Error z value Pr(>|z|)
> > (Intercept) -13.563 2.540 -5.339 9.33e-08 ***
> > riskItem -2.418 3.745 -0.646 0.518
> > ---
> > Signif. codes: 0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
> >
> > Correlation of Fixed Effects:
> > (Intr)
> > riskItem 0.064
> >
> > No effect but small dataset and test of principle works. However, when I
> > come to test the overall difference by items I hit
> >
> > 5 NA NA NA NA NA NA NA NA NA NA
> > NA NA NA NA NA NA
> > BITE23 BITE24 BITE25 BITE26 BITE27 BITE28 BITE29 BITE30 BITE31 BITE32
> > BITE33 Notas origRow Centro2 transfer one
> > 1 NA NA NA NA NA NA NA NA NA NA
> > NA <NA> 43 Urgell 0 1
> > 2 NA NA NA NA NA NA NA NA NA NA
> > NA <NA> 68 Urgell 0 1
> > 3 NA NA NA NA NA NA NA NA NA NA
> > NA <NA> 69 Urgell 0 1
> > 4 NA NA NA NA NA NA NA NA NA NA
> > NA <NA> 70 Urgell 0 1
> > 5 NA NA NA NA NA NA NA NA NA NA
> > NA <NA> 71 Urgell 0 1
> > Exclude quaireN ID.quaireN Nquaires ID.episodeN episodeN nEpisodes
> > Fecha.inicio Fecha.fin ExpCOREOM ExpSFA ExpSFB
> > 1 0 5 URGE0003-5 27 URGE0003-001 1 1
> > 2017-11-13 2019-02-12 NA NA NA
> > 2 0 3 URGE0004-3 36 URGE0004-001 1 2
> > 2017-11-21 2018-09-10 NA NA NA
> > 3 0 4 URGE0004-4 36 URGE0004-001 1 2
> > 2017-11-21 2018-09-10 NA NA NA
> > 4 0 5 URGE0004-5 36 URGE0004-001 1 2
> > 2017-11-21 2018-09-10 NA NA NA
> > 5 0 6 URGE0004-6 36 URGE0004-001 1 2
> > 2017-11-21 2018-09-10 NA NA NA
> > ExpYPCORE ExpEAT ExpBITE missingCOREOM missingSFA missingSFB missingYP
> > missingEAT missingBITE impossCOREOM impossSFA
> > 1 NA NA NA NA NA NA NA
> > NA NA NA NA
> > 2 NA NA NA NA NA NA NA
> > NA NA NA NA
> > 3 NA NA NA NA NA NA NA
> > NA NA NA NA
> > 4 NA NA NA NA NA NA NA
> > NA NA NA NA
> > 5 NA NA NA NA NA NA NA
> > NA NA NA NA
> > impossSFB impossYP impossEAT impossBITE
> > 1 NA NA NA NA
> > 2 NA NA NA NA
> > 3 NA NA NA NA
> > 4 NA NA NA NA
> > 5 NA NA NA NA
> > [ reached 'max' / getOption("max.print") -- omitted 672 rows ]
> > > datPrincip_URGE %>%
> > + arrange(Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + dim()
> > [1] 677 195
> > > datPrincip_URGE %>%
> > + arrange(Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + dim()
> > [1] 61 197
> > > dat <- bind_rows(datPrincip_ALCA,
> > + datPrincip_URGE)
> > > dat %>%
> > + arrange(Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + dim()
> > [1] 118 197
> > > dat <- bind_rows(datPrincip_ALCA,
> > + datPrincip_ARGE,
> > + datPrincip_AVEN,
> > + datPrincip_MOSC,
> > + datPrincip_SABA,
> > + datPrincip_SEVI,
> > + datPrincip_TARR,
> > + datPrincip_URGE)
> > > dat %>%
> > + arrange(Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + dim()
> > [1] 237 197
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + dim()
> > [1] 237 197
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + dim()
> > [1] 237 11
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + pivot_longer(cols = starts_with("YPCORE")) %>%
> > + dim()
> > [1] 2370 3
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + pivot_longer(cols = starts_with("YPCORE")) %>%
> > + mutate(missing = if_else(!is.na(value), 0, 1)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID name value missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > table(longDat$missing)
> > 0 1
> > 2356 14
> > > ?pivot_longer
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + pivot_longer(cols = starts_with("YPCORE"),
> > + names_to = "Item") %>%
> > + mutate(missing = if_else(!is.na(value), 0, 1)) -> longDat
> > > head(longD)
> > Error in head(longD) : object 'longD' not found
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID Item value missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + pivot_longer(cols = starts_with("YPCORE"),
> > + names_to = "Item",
> > + values_to = "Score") %>%
> > + mutate(missing = if_else(!is.na(Score), 0, 1)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > formul1 <- missing ~ Item + Item | ID
> > > fit <- glmer(formul1, family = binomial, data = longDat)
> > Error in glmer(formul1, family = binomial, data = longDat) :
> > could not find function "glmer"
> > > library(lme4)
> > Loading required package: Matrix
> >
> > Attaching package: ‘Matrix'
> >
> > The following objects are masked from ‘package:tidyr':
> >
> > expand, pack, unpack
> >
> > > formul1 <- missing ~ Item + Item | ID
> > > fit <- glmer(formul1, family = binomial, data = longDat)
> > boundary (singular) fit: see ?isSingular
> > Warning messages:
> > 1: In commonArgs(par, fn, control, environment()) :
> > maxfun < 10 * length(par)^2 is not recommended.
> > 2: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,
> > >
> > failure to converge in 10000 evaluations
> > > formul1 <- missing ~ Item | ID
> > > fit <- glmer(formul1, family = binomial, data = longDat)
> > boundary (singular) fit: see ?isSingular
> > Warning messages:
> > 1: In commonArgs(par, fn, control, environment()) :
> > maxfun < 10 * length(par)^2 is not recommended.
> > 2: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,
> > >
> > failure to converge in 10000 evaluations
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > ?glmer
> > > ?mcnemar.test
> > > table(is.na(dat$YPCORE01), is.na(dat$YPCORE02))
> > FALSE TRUE
> > FALSE 3573 4
> > TRUE 1 9499
> > > mcnemar.test(table(is.na(dat$YPCORE01), is.na(dat$YPCORE02)))
> > McNemar's Chi-squared test with continuity correction
> >
> > data: table(is.na(dat$YPCORE01), is.na(dat$YPCORE02))
> > McNemar's chi-squared = 0.8, df = 1, p-value = 0.3711
> >
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + pivot_longer(cols = c("YPCORE01", "YPCORE02")),
> > Error: unexpected ',' in:
> > " select(ID, starts_with("YP")) %>%
> > pivot_longer(cols = c("YPCORE01", "YPCORE02")),"
> > > names_to = "Item",
> > Error: unexpected ',' in " names_to = "Item","
> > > values_to = "Score") %>%
> > Error: unexpected ')' in " values_to = "Score")"
> > > mutate(missing = if_else(!is.na(Score), 0, 1)) -> longDat
> > Error in if_else(!is.na(Score), 0, 1) : object 'Score' not found
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YP")) %>%
> > + pivot_longer(cols = c("YPCORE01", "YPCORE02"),
> > + names_to = "Item",
> > + values_to = "Score") %>%
> > + mutate(missing = if_else(!is.na(Score), 0, 1)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 12
> > # Groups: ID [3]
> > ID YPCORE03 YPCORE04 YPCORE05 YPCORE06 YPCORE07 YPCORE08 YPCORE09
> > YPCORE10 Item Score missing
> > <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
> > <dbl> <chr> <dbl> <dbl>
> > 1 ALCA0001 0 0 0 1 0 0 0
> > 0 YPCORE01 1 0
> > 2 ALCA0001 0 0 0 1 0 0 0
> > 0 YPCORE02 0 0
> > 3 ALCA0004 0 0 0 1 0 0 0
> > 0 YPCORE01 1 0
> > 4 ALCA0004 0 0 0 1 0 0 0
> > 0 YPCORE02 0 0
> > 5 ALCA0007 1 0 0 1 0 2 1
> > 0 YPCORE01 2 0
> > 6 ALCA0007 1 0 0 1 0 2 1
> > 0 YPCORE02 1 0
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, c("YPCORE01", "YPCORE02")) %>%
> > + pivot_longer(cols = c("YPCORE01", "YPCORE02"),
> > + names_to = "Item",
> > + values_to = "Score") %>%
> > + mutate(missing = if_else(!is.na(Score), 0, 1)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [3]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0004 YPCORE01 1 0
> > 4 ALCA0004 YPCORE02 0 0
> > 5 ALCA0007 YPCORE01 2 0
> > 6 ALCA0007 YPCORE02 1 0
> > > fit <- glmer(formul1, family = binomial, data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > ?isSingular
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [3]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0004 YPCORE01 1 0
> > 4 ALCA0004 YPCORE02 0 0
> > 5 ALCA0007 YPCORE01 2 0
> > 6 ALCA0007 YPCORE02 1 0
> > > table(longDat$missing)
> > 0 1
> > 472 2
> > > library(lme4)
> > > formul1 <- missing ~ (Item | ID)
> > > fit <- glmer(formul1, family = binomial, data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > longDat$missRand <- rbinom(1, length(longDat), .5)
> > > head(longDat)
> > # A tibble: 6 x 5
> > # Groups: ID [3]
> > ID Item Score missing missRand
> > <chr> <chr> <dbl> <dbl> <int>
> > 1 ALCA0001 YPCORE01 1 0 2
> > 2 ALCA0001 YPCORE02 0 0 2
> > 3 ALCA0004 YPCORE01 1 0 2
> > 4 ALCA0004 YPCORE02 0 0 2
> > 5 ALCA0007 YPCORE01 2 0 2
> > 6 ALCA0007 YPCORE02 1 0 2
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > ?rbinom
> > > rbinom(1,1,.05)
> > [1] 0
> > > rbinom(1,1,.05)
> > [1] 0
> > > rbinom(1,1,.95)
> > [1] 1
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.95)
> > [1] 1
> > > rbinom(1,1,.5)
> > [1] 1
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 1
> > > rbinom(1,1,.5)
> > [1] 0
> > > rbinom(1,1,.5)
> > [1] 1
> > > rbinom(1,2,.5)
> > [1] 2
> > > rbinom(2,1,.5)
> > [1] 0 0
> > > rbinom(2,1,.5)
> > [1] 0 0
> > > rbinom(2,1,.5)
> > [1] 0 0
> > > rbinom(2,1,.5)
> > [1] 0 1
> > > rbinom(2,1,.5)
> > [1] 1 1
> > > rbinom(2,1,.5)
> > [1] 0 1
> > > longDat$missRand <- rbinom(length(longDat), 1, .5)
> > Error: Assigned data `rbinom(length(longDat), 1, 0.5)` must be compatible
> > with existing data.
> > x Existing data has 474 rows.
> > x Assigned data has 5 rows.
> > ℹ Only vectors of size 1 are recycled.
> > Run `rlang::last_error()` to see where the error occurred.
> > > table(longDat$missing)
> > 0 1
> > 472 2
> > > library(lme4)
> > > head(longDat)
> > # A tibble: 6 x 5
> > # Groups: ID [3]
> > ID Item Score missing missRand
> > <chr> <chr> <dbl> <dbl> <int>
> > 1 ALCA0001 YPCORE01 1 0 2
> > 2 ALCA0001 YPCORE02 0 0 2
> > 3 ALCA0004 YPCORE01 1 0 2
> > 4 ALCA0004 YPCORE02 0 0 2
> > 5 ALCA0007 YPCORE01 2 0 2
> > 6 ALCA0007 YPCORE02 1 0 2
> > > rbinom(2,1,.5)
> > [1] 1 0
> > > rbinom(2,1,.5)
> > [1] 1 1
> > > rbinom(2,1,.5)
> > [1] 1 1
> > > longDat$missRand <- NULL
> > > longDat$missRand <- rbinom(length(longDat), 1, .5)
> > Error: Assigned data `rbinom(length(longDat), 1, 0.5)` must be compatible
> > with existing data.
> > x Existing data has 474 rows.
> > x Assigned data has 4 rows.
> > ℹ Only vectors of size 1 are recycled.
> > Run `rlang::last_error()` to see where the error occurred.
> > > longDat$missRand <- rbinom(nrow(longDat), 1, .5)
> > > table(longDat$missRand)
> > 0 1
> > 240 234
> > > head(longDat)
> > # A tibble: 6 x 5
> > # Groups: ID [3]
> > ID Item Score missing missRand
> > <chr> <chr> <dbl> <dbl> <int>
> > 1 ALCA0001 YPCORE01 1 0 0
> > 2 ALCA0001 YPCORE02 0 0 1
> > 3 ALCA0004 YPCORE01 1 0 0
> > 4 ALCA0004 YPCORE02 0 0 1
> > 5 ALCA0007 YPCORE01 2 0 0
> > 6 ALCA0007 YPCORE02 1 0 0
> > > formul1 <- missRand ~ (Item | ID)
> > > fit <- glmer(formul1, family = binomial, data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > fit <- glmer(missRand ~ Item, family = binomial, data = longDat)
> > Error: No random effects terms specified in formula
> > > fit <- glmer(missRand ~ Item | ID, family = binomial, data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > fit <- glmer(factor(missRand) ~ Item | ID, family = binomial, data =
> > longDat)
> > boundary (singular) fit: see ?isSingular
> > > fit <- glmer(missRand ~ Item | ID, family = "logistic", data = longDat)
> > Error in get(family, mode = "function", envir = parent.frame(2)) :
> > object 'logistic' of mode 'function' was not found
> > > fit <- glmer(missRand ~ Item + (1 | ID), family = "binomial", data =
> > longDat)
> > boundary (singular) fit: see ?isSingular
> > > ?glmer
> > > longDat
> > # A tibble: 474 x 5
> > # Groups: ID [237]
> > ID Item Score missing missRand
> > <chr> <chr> <dbl> <dbl> <int>
> > 1 ALCA0001 YPCORE01 1 0 0
> > 2 ALCA0001 YPCORE02 0 0 1
> > 3 ALCA0004 YPCORE01 1 0 0
> > 4 ALCA0004 YPCORE02 0 0 1
> > 5 ALCA0007 YPCORE01 2 0 0
> > 6 ALCA0007 YPCORE02 1 0 0
> > 7 ALCA0012 YPCORE01 3 0 1
> > 8 ALCA0012 YPCORE02 1 0 1
> > 9 ALCA0013 YPCORE01 4 0 1
> > 10 ALCA0013 YPCORE02 3 0 1
> > # … with 464 more rows
> > > table(longDat$missRand, longDat$Item)
> > YPCORE01 YPCORE02
> > 0 125 115
> > 1 112 122
> > > fit <- glmer(missRand ~ factor(Item) + (1 | ID), family = "binomial",
> > data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > fit <- glmer(missRand ~ factor(Item) + (1 | factor(ID)), family =
> > "binomial", data = longDat)
> > Error: couldn't evaluate grouping factor factor(ID) within model frame:
> > try adding grouping factor to data frame explicitly if possible
> > > fit <- glmer(missRand ~ factor(Item) + (1 | ID), family = "binomial",
> > data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > rm(fit)
> > > fit <- glmer(missRand ~ factor(Item) + (1 | ID), family = "binomial",
> > data = longDat)
> > boundary (singular) fit: see ?isSingular
> > > fit
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: missRand ~ factor(Item) + (1 | ID)
> > Data: longDat
> > AIC BIC logLik deviance df.resid
> > 662.1833 674.6669 -328.0917 656.1833 471
> > Random effects:
> > Groups Name Std.Dev.
> > ID (Intercept) 1.943e-07
> > Number of obs: 474, groups: ID, 237
> > Fixed Effects:
> > (Intercept) factor(Item)YPCORE02
> > -0.1098 0.1689
> > convergence code 0; 0 optimizer warnings; 1 lme4 warnings
> > > dat %>%
> > + arrange(ID, Fecha.evaluacion) %>% # sort by evaluation date
> > + filter(Cuestionarios == "YP-CORE" &
> > + Codigo.datos == "RESPONDIO CORRECTAMENTE") %>%
> > + group_by(ID) %>%
> > + mutate(occasion = row_number(),
> > + first = if_else(occasion == 1, 1, 0)) %>% # to get first
> > ocassion only
> > + filter(first == 1) %>%
> > + select(ID, starts_with("YPCORE")) %>%
> > + pivot_longer(cols = starts_with("YPCORE"),
> > + names_to = "Item",
> > + values_to = "Score") %>%
> > + mutate(missing = if_else(!is.na(Score), 0, 1)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > formul1 <- missing ~ (Item | ID)
> > > fit <- glmer(missing ~ factor(Item) + (1 | ID), family = "binomial",
> > data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 1.78883 (tol = 0.002,
> > component 1)
> > > fit <- glmer((score == 1) ~ factor(Item) + (1 | ID), family =
> > "binomial", data = longDat)
> > Error in eval(predvars, data, env) : object 'score' not found
> > > fit <- glmer((Score == 1) ~ factor(Item) + (1 | ID), family =
> > "binomial", data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 0.00309555 (tol = 0.002,
> > component 1)
> > > fit <- glmer((Score == 1) ~ factor(Item) + (factor(Item) | ID), family =
> > "binomial", data = longDat)
> > Error: number of observations (=2356) < number of random effects (=2360)
> > for term (factor(Item) | ID); the random-effects parameters are probably
> > unidentifiable
> > > fit <- glmer((Score == 1) ~ factor(Item) + (Item | ID), family =
> > "binomial", data = longDat)
> > Error: number of observations (=2356) < number of random effects (=2360)
> > for term (Item | ID); the random-effects parameters are probably
> > unidentifiable
> > > fit <- glmer((Score == 1) ~ factor(Item) + (1 | ID), family =
> > "binomial", data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 0.00309555 (tol = 0.002,
> > component 1)
> > > head(longDat)
> > # A tibble: 6 x 4
> > # Groups: ID [1]
> > ID Item Score missing
> > <chr> <chr> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0
> > 2 ALCA0001 YPCORE02 0 0
> > 3 ALCA0001 YPCORE03 0 0
> > 4 ALCA0001 YPCORE04 0 0
> > 5 ALCA0001 YPCORE05 0 0
> > 6 ALCA0001 YPCORE06 1 0
> > > longDat$missRand <- rbinom(nrow(longDat), 1, .5)
> > > table(longDat$missRand)
> > 0 1
> > 1211 1159
> > > fit <- glmer(missRand ~ Item + (1 | ID), family = "binomial", data =
> > longDat)
> > boundary (singular) fit: see ?isSingular
> > > fit
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: missRand ~ Item + (1 | ID)
> > Data: longDat
> > AIC BIC logLik deviance df.resid
> > 3300.819 3364.296 -1639.410 3278.819 2359
> > Random effects:
> > Groups Name Std.Dev.
> > ID (Intercept) 2.663e-07
> > Number of obs: 2370, groups: ID, 237
> > Fixed Effects:
> > (Intercept) ItemYPCORE02 ItemYPCORE03 ItemYPCORE04 ItemYPCORE05
> > ItemYPCORE06 ItemYPCORE07 ItemYPCORE08
> > -5.909e-02 1.013e-01 -1.690e-02 1.182e-01 6.753e-02
> > -1.186e-01 -1.046e-15 -1.186e-01
> > ItemYPCORE09 ItemYPCORE10
> > 1.858e-01 -6.766e-02
> > convergence code 0; 0 optimizer warnings; 1 lme4 warnings
> > > anova(fit)
> > Analysis of Variance Table
> > npar Sum Sq Mean Sq F value
> > Item 9 5.5453 0.61615 0.6161
> > > table(longDat$Item, as.numeric(longDat$Item))
> > < table of extent 10 x 0 >
> > Warning message:
> > In table(longDat$Item, as.numeric(longDat$Item)) :
> > NAs introduced by coercion
> > > dim(longDat)
> > [1] 2370 5
> > > head(longDat)
> > # A tibble: 6 x 5
> > # Groups: ID [1]
> > ID Item Score missing missRand
> > <chr> <chr> <dbl> <dbl> <int>
> > 1 ALCA0001 YPCORE01 1 0 1
> > 2 ALCA0001 YPCORE02 0 0 1
> > 3 ALCA0001 YPCORE03 0 0 1
> > 4 ALCA0001 YPCORE04 0 0 1
> > 5 ALCA0001 YPCORE05 0 0 1
> > 6 ALCA0001 YPCORE06 1 0 0
> > > table(longDat$Item) #, as.numeric(longDat$Item))
> > YPCORE01 YPCORE02 YPCORE03 YPCORE04 YPCORE05 YPCORE06 YPCORE07 YPCORE08
> > YPCORE09 YPCORE10
> > 237 237 237 237 237 237 237 237
> > 237 237
> > > table(as.numeric(longDat$Item)) #, as.numeric(longDat$Item))
> > < table of extent 0 >
> > Warning message:
> > In table(as.numeric(longDat$Item)) : NAs introduced by coercion
> > > table(longDat$Item, as.numeric(factor(longDat$Item)))
> > 1 2 3 4 5 6 7 8 9 10
> > YPCORE01 237 0 0 0 0 0 0 0 0 0
> > YPCORE02 0 237 0 0 0 0 0 0 0 0
> > YPCORE03 0 0 237 0 0 0 0 0 0 0
> > YPCORE04 0 0 0 237 0 0 0 0 0 0
> > YPCORE05 0 0 0 0 237 0 0 0 0 0
> > YPCORE06 0 0 0 0 0 237 0 0 0 0
> > YPCORE07 0 0 0 0 0 0 237 0 0 0
> > YPCORE08 0 0 0 0 0 0 0 237 0 0
> > YPCORE09 0 0 0 0 0 0 0 0 237 0
> > YPCORE10 0 0 0 0 0 0 0 0 0 237
> > > longDat %>%
> > + mutate(missRand = rbinom(1, 1, .5))
> > # A tibble: 2,370 x 5
> > # Groups: ID [237]
> > ID Item Score missing missRand
> > <chr> <chr> <dbl> <dbl> <int>
> > 1 ALCA0001 YPCORE01 1 0 1
> > 2 ALCA0001 YPCORE02 0 0 1
> > 3 ALCA0001 YPCORE03 0 0 1
> > 4 ALCA0001 YPCORE04 0 0 1
> > 5 ALCA0001 YPCORE05 0 0 1
> > 6 ALCA0001 YPCORE06 1 0 1
> > 7 ALCA0001 YPCORE07 0 0 1
> > 8 ALCA0001 YPCORE08 0 0 1
> > 9 ALCA0001 YPCORE09 0 0 1
> > 10 ALCA0001 YPCORE10 0 0 1
> > # … with 2,360 more rows
> > > table(longDat$missRand)
> > 0 1
> > 1211 1159
> > > longDat %>%
> > + mutate(missRand = rbinom(1, 1, .5),
> > + missReg = rbinom(1, 1, as.numeric(factor(Item))/11))
> > # A tibble: 2,370 x 6
> > # Groups: ID [237]
> > ID Item Score missing missRand missReg
> > <chr> <chr> <dbl> <dbl> <int> <int>
> > 1 ALCA0001 YPCORE01 1 0 1 0
> > 2 ALCA0001 YPCORE02 0 0 1 0
> > 3 ALCA0001 YPCORE03 0 0 1 0
> > 4 ALCA0001 YPCORE04 0 0 1 0
> > 5 ALCA0001 YPCORE05 0 0 1 0
> > 6 ALCA0001 YPCORE06 1 0 1 0
> > 7 ALCA0001 YPCORE07 0 0 1 0
> > 8 ALCA0001 YPCORE08 0 0 1 0
> > 9 ALCA0001 YPCORE09 0 0 1 0
> > 10 ALCA0001 YPCORE10 0 0 1 0
> > # … with 2,360 more rows
> > > table(longDat$missReg)
> > < table of extent 0 >
> > Warning message:
> > Unknown or uninitialised column: `missReg`.
> > > longDat %>%
> > + mutate(missRand = rbinom(1, 1, .5),
> > + missReg = rbinom(1, 1, as.numeric(factor(Item))/11)) -> longDat
> > > table(longDat$missReg)
> > 0 1
> > 2120 250
> > > fit <- glmer(missReg ~ Item + (1 | ID), family = "binomial", data =
> > longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 3.81902 (tol = 0.002,
> > component 1)
> > > longDat %>%
> > + mutate(missRand = rbinom(1, 1, .5),
> > + missReg = rbinom(1, 1, as.numeric(factor(Item))/11)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 6
> > # Groups: ID [1]
> > ID Item Score missing missRand missReg
> > <chr> <chr> <dbl> <dbl> <int> <int>
> > 1 ALCA0001 YPCORE01 1 0 1 1
> > 2 ALCA0001 YPCORE02 0 0 1 1
> > 3 ALCA0001 YPCORE03 0 0 1 1
> > 4 ALCA0001 YPCORE04 0 0 1 1
> > 5 ALCA0001 YPCORE05 0 0 1 1
> > 6 ALCA0001 YPCORE06 1 0 1 1
> > > table(longDat$missRand)
> > 0 1
> > 1140 1230
> > > table(longDat$missReg)
> > 0 1
> > 2160 210
> > > library(lme4)
> > > fitRand <- glmer(missRand ~ Item + (1 | ID), family = "binomial", data =
> > longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 5.46528 (tol = 0.002,
> > component 1)
> > > anova(fitRand)
> > Analysis of Variance Table
> > npar Sum Sq Mean Sq F value
> > Item 9 0.050372 0.0055969 0.0056
> > > length(unique(longDat$ID))
> > [1] 237
> > > rm(fitRand)
> > > anova(fitRand)
> > Error in anova(fitRand) : object 'fitRand' not found
> > > fitRand <- glmer(missRand ~ Item + (1 | ID), family = "binomial", data =
> > longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 5.46528 (tol = 0.002,
> > component 1)
> > > fitRand <- glmer(missRand ~ Item + (1 | ID), family = binomial, data =
> > longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 5.46528 (tol = 0.002,
> > component 1)
> > > fitRand <- glmer(missRand ~ Item + (1 | ID), family = binomial(link =
> > "logit"), data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 5.46528 (tol = 0.002,
> > component 1)
> > > ?glmer
> > > data("cbpp")
> > > dim(cbpp)
> > [1] 56 4
> > > head(cbpp)
> > herd incidence size period
> > 1 1 2 14 1
> > 2 1 3 12 2
> > 3 1 4 9 3
> > 4 1 0 5 4
> > 5 2 3 22 1
> > 6 2 1 18 2
> > > fitRand <- glmer(cbind(missRand, 1) ~ Item + (1 | ID), family =
> > binomial(link = "logit"), data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 0.00550177 (tol = 0.002,
> > component 1)
> > > fitReg <- glmer(cbind(missReg, 1) ~ Item + (1 | ID), family =
> > "binomial", data = longDat)
> > Warning messages:
> > 1: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,
> > >
> > failure to converge in 10000 evaluations
> > 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 1.47499 (tol = 0.002,
> > component 1)
> > > fitReal <- glmer(cbind(missing, 1) ~ Item + (1 | ID), family =
> > "binomial", data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 1.29013 (tol = 0.002,
> > component 1)
> > > data(mlmdata)
> > Warning message:
> > In data(mlmdata) : data set ‘mlmdata' not found
> > > library(haven)
> > > mlmdata <- read_dta("
> > https://stats.idre.ucla.edu/stat/examples/imm/imm10.dta")
> > > head(mlmdata)
> > # A tibble: 6 x 19
> > schid stuid ses meanses homework white parented public ratio percmin
> > math sex race sctype cstr scsize urban
> > <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
> > <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
> > 1 7472 3 -0.130 -0.483 1 1 2 1 19 0
> > 48 2 4 1 2 3 2
> > 2 7472 8 -0.390 -0.483 0 1 2 1 19 0
> > 48 1 4 1 2 3 2
> > 3 7472 13 -0.800 -0.483 0 1 2 1 19 0
> > 53 1 4 1 2 3 2
> > 4 7472 17 -0.720 -0.483 1 1 2 1 19 0
> > 42 1 4 1 2 3 2
> > 5 7472 27 -0.740 -0.483 2 1 2 1 19 0
> > 43 2 4 1 2 3 2
> > 6 7472 28 -0.580 -0.483 1 1 2 1 19 0
> > 57 2 4 1 2 3 2
> > # … with 2 more variables: region <dbl>, schnum <dbl>
> > > model <- glmer(white ~ homework + (1 + homework | schid), data=mlmdata,
> > + family=binomial(link="logit"))
> > boundary (singular) fit: see ?isSingular
> > > summary(model)
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: white ~ homework + (1 + homework | schid)
> > Data: mlmdata
> >
> > AIC BIC logLik deviance df.resid
> > 182.4 200.2 -86.2 172.4 255
> >
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -4.3373 -0.1184 0.1112 0.3421 3.8801
> >
> > Random effects:
> > Groups Name Variance Std.Dev. Corr
> > schid (Intercept) 16.28745 4.0358
> > homework 0.04678 0.2163 -1.00
> > Number of obs: 260, groups: schid, 10
> >
> > Fixed effects:
> > Estimate Std. Error z value Pr(>|z|)
> > (Intercept) 1.67365 1.47324 1.136 0.256
> > homework 0.04508 0.18433 0.245 0.807
> >
> > Correlation of Fixed Effects:
> > (Intr)
> > homework -0.590
> > convergence code: 0
> > boundary (singular) fit: see ?isSingular
> >
> > > dim(mlmdata)
> > [1] 260 19
> > > head(mlmdata)
> > # A tibble: 6 x 19
> > schid stuid ses meanses homework white parented public ratio percmin
> > math sex race sctype cstr scsize urban
> > <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
> > <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
> > 1 7472 3 -0.130 -0.483 1 1 2 1 19 0
> > 48 2 4 1 2 3 2
> > 2 7472 8 -0.390 -0.483 0 1 2 1 19 0
> > 48 1 4 1 2 3 2
> > 3 7472 13 -0.800 -0.483 0 1 2 1 19 0
> > 53 1 4 1 2 3 2
> > 4 7472 17 -0.720 -0.483 1 1 2 1 19 0
> > 42 1 4 1 2 3 2
> > 5 7472 27 -0.740 -0.483 2 1 2 1 19 0
> > 43 2 4 1 2 3 2
> > 6 7472 28 -0.580 -0.483 1 1 2 1 19 0
> > 57 2 4 1 2 3 2
> > # … with 2 more variables: region <dbl>, schnum <dbl>
> > > length(unique(mlmdata$schid))
> > [1] 10
> > > table(white)
> > Error in table(white) : object 'white' not found
> > > table(mlmdata$white)
> > 0 1
> > 71 189
> > > fitRand <- glmer(missRand ~ Item + (1 | ID), family = binomial(link =
> > "logit"), data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 5.46528 (tol = 0.002,
> > component 1)
> > > ,
> > Error: unexpected ',' in ","
> > > longDat %>%
> > + mutate(missRand = rbinom(1, 1, .5),
> > + missReg = rbinom(1, 1, as.numeric(factor(Item))/11),
> > + riskItem = if_else(Item == "YPCORE04", 1, 0),
> > + posItem = if_else(Item %in% c("YPCORE03", "YPCORE5",
> > "YPCORE10"), 1, 0)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 8
> > # Groups: ID [1]
> > ID Item Score missing missRand missReg riskItem posItem
> > <chr> <chr> <dbl> <dbl> <int> <int> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0 0 0 0 0
> > 2 ALCA0001 YPCORE02 0 0 0 0 0 0
> > 3 ALCA0001 YPCORE03 0 0 0 0 0 1
> > 4 ALCA0001 YPCORE04 0 0 0 0 1 0
> > 5 ALCA0001 YPCORE05 0 0 0 0 0 0
> > 6 ALCA0001 YPCORE06 1 0 0 0 0 0
> > > table(longDat$riskItem, longDat$posItem)
> > 0 1
> > 0 1659 474
> > 1 237 0
> > > table(longDat$riskItem, longDat$Item)
> > YPCORE01 YPCORE02 YPCORE03 YPCORE04 YPCORE05 YPCORE06 YPCORE07
> > YPCORE08 YPCORE09 YPCORE10
> > 0 237 237 237 0 237 237 237
> > 237 237 237
> > 1 0 0 0 237 0 0 0
> > 0 0 0
> > > table(longDat$posItem, longDat$Item)
> > YPCORE01 YPCORE02 YPCORE03 YPCORE04 YPCORE05 YPCORE06 YPCORE07
> > YPCORE08 YPCORE09 YPCORE10
> > 0 237 237 0 237 237 237 237
> > 237 237 0
> > 1 0 0 237 0 0 0 0
> > 0 0 237
> > > longDat %>%
> > + mutate(missRand = rbinom(1, 1, .5),
> > + missReg = rbinom(1, 1, as.numeric(factor(Item))/11),
> > + riskItem = if_else(Item == "YPCORE04", 1, 0),
> > + posItem = if_else(Item %in% c("YPCORE03", "YPCORE05",
> > "YPCORE10"), 1, 0)) -> longDat
> > > head(longDat)
> > # A tibble: 6 x 8
> > # Groups: ID [1]
> > ID Item Score missing missRand missReg riskItem posItem
> > <chr> <chr> <dbl> <dbl> <int> <int> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0 0 0 0 0
> > 2 ALCA0001 YPCORE02 0 0 0 0 0 0
> > 3 ALCA0001 YPCORE03 0 0 0 0 0 1
> > 4 ALCA0001 YPCORE04 0 0 0 0 1 0
> > 5 ALCA0001 YPCORE05 0 0 0 0 0 1
> > 6 ALCA0001 YPCORE06 1 0 0 0 0 0
> > > table(longDat$posItem, longDat$Item)
> > YPCORE01 YPCORE02 YPCORE03 YPCORE04 YPCORE05 YPCORE06 YPCORE07
> > YPCORE08 YPCORE09 YPCORE10
> > 0 237 237 0 237 0 237 237
> > 237 237 0
> > 1 0 0 237 0 237 0 0
> > 0 0 237
> > > fitRisk <- glmer(missing ~ riskItem + (1 | ID), family = binomial, data
> > = longDat)
> > > summary(fitRisk)
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: missing ~ riskItem + (1 | ID)
> > Data: longDat
> >
> > AIC BIC logLik deviance df.resid
> > 51.7 69.0 -22.9 45.7 2367
> >
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -0.52020 -0.00113 -0.00113 -0.00113 2.93829
> >
> > Random effects:
> > Groups Name Variance Std.Dev.
> > ID (Intercept) 205.3 14.33
> > Number of obs: 2370, groups: ID, 237
> >
> > Fixed effects:
> > Estimate Std. Error z value Pr(>|z|)
> > (Intercept) -13.563 2.540 -5.339 9.33e-08 ***
> > riskItem -2.418 3.745 -0.646 0.518
> > ---
> > Signif. codes: 0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
> >
> > Correlation of Fixed Effects:
> > (Intr)
> > riskItem 0.064
> > > fitRisk <- glmer(missing ~ posItem + (1 | ID), family = binomial, data =
> > longDat)
> > > fitRisk <- glmer(missing ~ riskItem + (1 | ID), family = binomial, data
> > = longDat)
> > > summary(fitRisk)
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: missing ~ riskItem + (1 | ID)
> > Data: longDat
> >
> > AIC BIC logLik deviance df.resid
> > 51.7 69.0 -22.9 45.7 2367
> >
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -0.52020 -0.00113 -0.00113 -0.00113 2.93829
> >
> > Random effects:
> > Groups Name Variance Std.Dev.
> > ID (Intercept) 205.3 14.33
> > Number of obs: 2370, groups: ID, 237
> >
> > Fixed effects:
> > Estimate Std. Error z value Pr(>|z|)
> > (Intercept) -13.563 2.540 -5.339 9.33e-08 ***
> > riskItem -2.418 3.745 -0.646 0.518
> > ---
> > Signif. codes: 0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
> >
> > Correlation of Fixed Effects:
> > (Intr)
> > riskItem 0.064
> > > fitPos <- glmer(missing ~ posItem + (1 | ID), family = binomial, data =
> > longDat)
> > > summary(fitPos)
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: missing ~ posItem + (1 | ID)
> > Data: longDat
> >
> > AIC BIC logLik deviance df.resid
> > 52.4 69.7 -23.2 46.4 2367
> >
> > Scaled residuals:
> > Min 1Q Median 3Q Max
> > -0.5111 -0.0012 -0.0012 -0.0010 3.4498
> >
> > Random effects:
> > Groups Name Variance Std.Dev.
> > ID (Intercept) 197.5 14.05
> > Number of obs: 2370, groups: ID, 237
> >
> > Fixed effects:
> > Estimate Std. Error z value Pr(>|z|)
> > (Intercept) -13.5152 2.5280 -5.346 8.98e-08 ***
> > posItem -0.2953 1.2396 -0.238 0.812
> > ---
> > Signif. codes: 0 ‘***' 0.001 ‘**' 0.01 ‘*' 0.05 ‘.' 0.1 ‘ ' 1
> >
> > Correlation of Fixed Effects:
> > (Intr)
> > posItem -0.110
> > > head(longDat,15)
> > # A tibble: 15 x 8
> > # Groups: ID [2]
> > ID Item Score missing missRand missReg riskItem posItem
> > <chr> <chr> <dbl> <dbl> <int> <int> <dbl> <dbl>
> > 1 ALCA0001 YPCORE01 1 0 0 0 0 0
> > 2 ALCA0001 YPCORE02 0 0 0 0 0 0
> > 3 ALCA0001 YPCORE03 0 0 0 0 0 1
> > 4 ALCA0001 YPCORE04 0 0 0 0 1 0
> > 5 ALCA0001 YPCORE05 0 0 0 0 0 1
> > 6 ALCA0001 YPCORE06 1 0 0 0 0 0
> > 7 ALCA0001 YPCORE07 0 0 0 0 0 0
> > 8 ALCA0001 YPCORE08 0 0 0 0 0 0
> > 9 ALCA0001 YPCORE09 0 0 0 0 0 0
> > 10 ALCA0001 YPCORE10 0 0 0 0 0 1
> > 11 ALCA0004 YPCORE01 1 0 0 0 0 0
> > 12 ALCA0004 YPCORE02 0 0 0 0 0 0
> > 13 ALCA0004 YPCORE03 0 0 0 0 0 1
> > 14 ALCA0004 YPCORE04 0 0 0 0 1 0
> > 15 ALCA0004 YPCORE05 0 0 0 0 0 1
> > >
> > > rm(fitRand)
> > > rm(fitReg)
> > > rm(fitReal)
> > > fitRand <- glmer(cbind(missRand, 1) ~ Item + (1 | ID), family =
> > binomial(link = "logit"), data = longDat)
> > Warning message:
> > In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 0.0501701 (tol = 0.002,
> > component 1)
> > > fitRand
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: cbind(missRand, 1) ~ Item + (1 | ID)
> > Data: longDat
> > AIC BIC logLik deviance df.resid
> > 2215.017 2278.494 -1096.509 2193.017 2359
> > Random effects:
> > Groups Name Std.Dev.
> > ID (Intercept) 2.509
> > Number of obs: 2370, groups: ID, 237
> > Fixed Effects:
> > (Intercept) ItemYPCORE02 ItemYPCORE03 ItemYPCORE04 ItemYPCORE05
> > ItemYPCORE06 ItemYPCORE07 ItemYPCORE08
> > -2.512456 0.004360 0.007166 0.001288 0.004078
> > 0.002322 0.005059 0.006118
> > ItemYPCORE09 ItemYPCORE10
> > 0.006454 0.005057
> > convergence code 0; 0 optimizer warnings; 1 lme4 warnings
> >
> > So failed to converge but does give a set of estimates. I think I can use
> > the
> >
> >
> > > fitReg <- glmer(cbind(missReg, 1) ~ Item + (1 | ID), family =
> > "binomial", data = longDat)
> > Warning messages:
> > 1: In (function (fn, par, lower = rep.int(-Inf, n), upper = rep.int(Inf,
> > >
> > failure to converge in 10000 evaluations
> > 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
> > Model failed to converge with max|grad| = 1.66664 (tol = 0.002,
> > component 1)
> > > fitReg
> > Generalized linear mixed model fit by maximum likelihood (Laplace
> > Approximation) ['glmerMod']
> > Family: binomial ( logit )
> > Formula: cbind(missReg, 1) ~ Item + (1 | ID)
> > Data: longDat
> > AIC BIC logLik deviance df.resid
> > 605.3404 668.8175 -291.6702 583.3404 2359
> > Random effects:
> > Groups Name Std.Dev.
> > ID (Intercept) 6.775
> > Number of obs: 2370, groups: ID, 237
> > Fixed Effects:
> > (Intercept) ItemYPCORE02 ItemYPCORE03 ItemYPCORE04 ItemYPCORE05
> > ItemYPCORE06 ItemYPCORE07 ItemYPCORE08
> > -8.6607 -0.3793 -0.3159 -0.3591 -0.4608
> > -0.4453 -0.4165 -0.6597
> > ItemYPCORE09 ItemYPCORE10
> > -0.4304 -0.3535
> > convergence code 0; 1 optimizer warnings; 1 lme4 warnings
> >
> > I am starting to realise the depths of my ignorance about all this. It's
> > clear to me that having only one observation per cell
> > because of the nesting of items within participants and only having one
> > completion per participant is causing the convergence
> > issues (which makes sense though I have a sense that there might be ways
> > to extract estimates despite this: this is beyond me!)
> >
> > I'm puzzled that I get slightly different results for the risk analysis if
> > I use the "cbind(missing, 1) ~ " syntax in the formula
> > from those I get using just "missing ~ " and different max|grad values
> > from the nonconvergence message depending on that choice.
> > I suspect that's a red herring.
> >
> > I have seen many comments here recently about non-convergence and ways to
> > overcome it by specifying different methods of
> > approximation (if I understood that properly) but they generally seemed to
> > be for much more complex models than mine and
> > probably weren't about this "cell size = 1" issue. I have searched for
> > answer or illumination for some hours and perhaps
> > I'm asking the wrong questions but I'm stuck. (I think it's very unlikely
> > that I have access to trained statisticians by
> > virtue of my honorary or paid academic positions!)
> >
> > So my questions:
> > 1) _is_ there a way to estimate such a model, i.e. estimating an effect
> > across all ten items with only one completion per participant?
> > 2) can someone suggest good reading matter for a dogged non-statistician
> > who can handle a lot of continuous variable issues up to
> > some real complexity on his own but has never been as comfortable with
> > counts beyond the basic old NHST, single level approaches? I
> > have an old copy of Pinheiro and Bates but confess that despite trying
> > many times, the older it and I get, the less I'm able to cope
> > with it. The algebra is just beyond me.
> >
> > TIA ... and best wishes to all for physical safety and psychological
> > resilience in these times,
> >
> > Chris
> >
> > --
> > Small contribution in our coronavirus rigours:
> >
> > https://www.coresystemtrust.org.uk/home/free-options-to-replace-paper-core-forms-during-the-coronavirus-pandemic/
> >
> > Chris Evans <chris@psyctc.org> Visiting Professor, University of
> > Sheffield <chris.evans@sheffield.ac.uk>
> > I do some consultation work for the University of Roehampton <
> > chris.evans@roehampton.ac.uk> and other places
> > but <chris@psyctc.org> remains my main Email address. I have a work web
> > site at:
> > https://www.psyctc.org/psyctc/
> > and a site I manage for CORE and CORE system trust at:
> > http://www.coresystemtrust.org.uk/
> > I have "semigrated" to France, see:
> > https://www.psyctc.org/pelerinage2016/semigrating-to-france/
> >
> > https://www.psyctc.org/pelerinage2016/register-to-get-updates-from-pelerinage2016/
> >
> > If you want an Emeeting, I am trying to keep them to Thursdays and my
> > diary is at:
> > https://www.psyctc.org/pelerinage2016/ceworkdiary/
> > Beware: French time, generally an hour ahead of UK.
> >
> > _______________________________________________
> > R-sig-mixed-models@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models@r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
_______________________________________________
R-sig-mixed-models@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[prev in list] [next in list] [prev in thread] [next in thread]
Configure |
About |
News |
Add a list |
Sponsored by KoreLogic