[R-sig-ME] Struggling (probably problem in wetware not software) with multilevel logistic regression
Phillip Alday
ph||||p@@|d@y @end|ng |rom mp|@n|
Mon May 4 01:11:27 CEST 2020
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 using 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 using psyctc.org> Visiting Professor, University of
>> Sheffield <chris.evans using sheffield.ac.uk>
>> I do some consultation work for the University of Roehampton <
>> chris.evans using roehampton.ac.uk> and other places
>> but <chris using 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 using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
More information about the R-sig-mixed-models
mailing list