1
$\begingroup$

I have two-time point data, where I want to measure "change" in a dependent variable (si_ch_r) predicted by change in three independent variables and their two- and three-way-interaction effects (pb_ch, tb_ch, bhs_ch). I have used raw-change / difference-scores to model "change".

I want to use a mixed-effect model so that missing values can be handled by maximum-likelihood estimation (my assumption is that this is a feasible approach). In order to get the model to run, I have specified the participant identifier in the random intercept term ((1 | ID)). My questions are:

  1. Is it feasible to specify ID as a random intercept even though each participant only has one observation in the model (due to the use of raw-change scores)?
  2. In the model summary, does Number of obs: 20, groups: ID, 20 mean that this is the number of observations with complete data, or is it actually dropping participants from the analysis (sample N = 30)?

I have included some basic descriptive statistics, my code, and the output of the fitted model below.

Please let me know if there is anything I have not provided in order to make answering my questions easier.

Thanks in advance.

library(psych)
library(glmmTMB)
library(tidyverse)

lpa_t0t1_outc_w <-
structure(list(si_ch_r = c(0.896473364001916, NA, NA, 0.911560735067541,
0.60376483779863, 0.219834563805069, 0.866894166636438, 0.63412384852177,
NA, NA, -0.696684917063051, -1.69539771027214, -0.0125334695080693,
NA, 0.415193850778427, NA, -0.98220269533347, NA, 0.615840188747972,
0.339809491013167, 0.240426031142308, NA, -1.16011988299752,
-0.240426031142308, -0.768820293458062, -0.91918273514682, -1.71688601843104,
-1.06693763219277, 0.874217164866484, -0.966088297132373), pb_ch = c(12,
5, -13, 14, 3, 18, 10, -1, 4, 2, 12, -2, -6, NA, 4, 8, NA, 0,
0, 9, 15, 6, -3, 7, 1, 9, 4, 5, 22, 13), tb_ch = c(11, 2, -2,
21, 6, 19, 15, 4, 3, 6, 14, -7, 7, NA, 9, 15, NA, -9, 4, 5, 11,
6, 5, -2, 3, 19, 1, 6, -7.75, 29), bhs_ch = c(10, -1, -2, 6,
10, 9, 8, 3, 1, 5, 3.75, 0, 0, NA, 1, 4, NA, 3, 10, 10, 8, 0,
NA, -1, -1, 5, 5, 6, 8, 11), ID = structure(c(351L, 429L, 483L,
389L, 201L, 210L, 83L, 279L, 127L, 357L, 336L, 139L, 374L, 65L,
183L, 334L, 97L, 467L, 425L, 245L, 456L, 41L, 464L, 217L, 113L,
430L, 56L, 492L, 345L, 312L), levels = c("2", "4", "6", "8",
"9", "10", "11", "12", "13", "14", "16", "17", "18", "19", "20",
"21", "22", "23", "24", "25", "26", "27", "30", "31", "33", "34",
"35", "36", "40", "42", "43", "44", "45", "48", "50", "51", "52",
"53", "54", "56", "57", "59", "60", "64", "65", "66", "67", "72",
"73", "74", "75", "76", "77", "78", "79", "80", "84", "85", "86",
"89", "90", "93", "94", "95", "96", "97", "99", "101", "102",
"103", "104", "105", "106", "107", "109", "110", "111", "112",
"113", "114", "115", "116", "117", "118", "119", "120", "121",
"123", "124", "125", "126", "127", "128", "129", "130", "131",
"132", "134", "135", "136", "137", "138", "139", "140", "141",
"142", "143", "144", "146", "148", "149", "150", "151", "152",
"153", "155", "156", "157", "159", "162", "163", "164", "165",
"166", "167", "168", "170", "173", "174", "176", "178", "182",
"183", "185", "186", "187", "189", "191", "192", "196", "197",
"198", "199", "201", "202", "203", "204", "205", "206", "207",
"208", "209", "212", "213", "214", "215", "216", "217", "218",
"222", "223", "224", "225", "226", "228", "231", "234", "235",
"236", "237", "239", "240", "241", "242", "244", "246", "247",
"250", "252", "253", "254", "255", "257", "258", "259", "260",
"261", "262", "264", "266", "267", "268", "269", "270", "271",
"272", "274", "275", "276", "277", "278", "279", "280", "282",
"284", "286", "287", "288", "293", "294", "295", "297", "298",
"299", "300", "301", "302", "303", "304", "306", "307", "308",
"309", "310", "311", "314", "317", "318", "321", "323", "327",
"331", "332", "333", "337", "338", "341", "342", "343", "344",
"345", "347", "349", "350", "351", "353", "354", "356", "357",
"358", "359", "362", "364", "366", "367", "371", "372", "377",
"378", "379", "380", "381", "382", "383", "385", "386", "387",
"393", "394", "395", "398", "400", "402", "404", "407", "410",
"411", "412", "413", "414", "415", "416", "417", "418", "419",
"420", "421", "422", "423", "427", "428", "429", "432", "433",
"434", "436", "437", "438", "440", "441", "444", "446", "448",
"449", "451", "453", "454", "455", "456", "457", "458", "459",
"461", "462", "463", "464", "465", "467", "468", "470", "473",
"474", "476", "477", "478", "479", "480", "481", "482", "483",
"484", "486", "487", "490", "491", "497", "499", "503", "504",
"505", "506", "507", "508", "509", "510", "513", "514", "518",
"522", "523", "524", "525", "526", "529", "530", "531", "532",
"533", "534", "535", "536", "537", "539", "541", "542", "543",
"545", "546", "547", "548", "550", "551", "553", "554", "556",
"557", "559", "560", "561", "562", "564", "565", "566", "567",
"568", "570", "572", "574", "575", "577", "578", "579", "580",
"582", "583", "584", "585", "586", "587", "588", "591", "592",
"593", "595", "596", "597", "601", "602", "603", "604", "605",
"606", "611", "613", "615", "616", "617", "618", "620", "621",
"622", "623", "624", "625", "626", "628", "632", "633", "634",
"636", "637", "640", "644", "645", "647", "648", "649", "650",
"651", "652", "653", "655", "656", "657", "659", "660", "661",
"662", "665", "666", "667", "668", "669", "670", "673", "675",
"676", "677", "678", "681", "682", "683", "684", "686", "690",
"691", "695", "696", "698", "700", "701", "702", "704", "705",
"710", "714", "715", "716", "717", "719", "722", "723", "724",
"725", "726", "728", "729", "730", "732", "733", "734", "736",
"738", "739", "741", "743", "745", "746", "749", "750"), class = "factor")), row.names = c(NA,
-30L), class = "data.frame")

lpa_t0t1_outc_w %>%
  select(si_ch_r, pb_ch, tb_ch, bhs_ch, ID) %>%
  describe()
#>         vars  n   mean     sd median trimmed    mad    min    max  range  skew
#> si_ch_r    1 22  -0.16   0.88    0.1   -0.11   1.16  -1.72   0.91   2.63 -0.29
#> pb_ch      2 28   5.64   7.54    5.0    5.71   7.41 -13.00  22.00  35.00 -0.10
#> tb_ch      3 28   6.90   8.89    6.0    6.67   7.41  -9.00  29.00  38.00  0.35
#> bhs_ch     4 27   4.51   4.16    5.0    4.51   5.93  -2.00  11.00  13.00  0.01
#> ID*        5 30 283.33 145.87  323.0  287.33 174.21  41.00 492.00 451.00 -0.20
#>         kurtosis    se
#> si_ch_r    -1.44  0.19
#> pb_ch      -0.05  1.43
#> tb_ch      -0.17  1.68
#> bhs_ch     -1.46  0.80
#> ID*        -1.42 26.63

lpa_t0t1_outc_w %>%
  select(ID, si_ch_r, pb_ch, tb_ch, bhs_ch) %>%
  pivot_longer(-ID) %>%
  group_by(ID) %>%
  mutate(miss = sum(value)) %>%
  slice(1) %>%
  summarise(is_na = sum(is.na(miss))) %>%
  filter(is_na %in% 1) %>%
  summarise(miss_any = sum(is_na)) %>%
  mutate(obs_all = 30 - miss_any)
#> # A tibble: 1 × 2
#>   miss_any obs_all
#>      <int>   <dbl>
#> 1       10      20

full <-
  glmmTMB(
    si_ch_r ~ pb_ch * tb_ch * bhs_ch + (1 | ID),
    na.action = na.exclude,
    data = lpa_t0t1_outc_w,
    family = gaussian,
    ziformula = ~0,
    REML = FALSE,
    control = glmmTMBControl(
      optCtrl = list(iter.max = 30000, eval.max = 40000),
      profile = TRUE
    )
  )

summary(full)
#>  Family: gaussian  ( identity )
#> Formula:          si_ch_r ~ pb_ch * tb_ch * bhs_ch + (1 | ID)
#> Data: lpa_t0t1_outc_w
#> 
#>      AIC      BIC   logLik deviance df.resid 
#>     61.6     71.5    -20.8     41.6       10 
#> 
#> Random effects:
#> 
#> Conditional model:
#>  Groups   Name        Variance Std.Dev.
#>  ID       (Intercept) 0.09371  0.3061  
#>  Residual             0.37466  0.6121  
#> Number of obs: 20, groups:  ID, 20
#> 
#> Dispersion estimate for gaussian family (sigma^2): 0.375 
#> 
#> Conditional model:
#>                      Estimate Std. Error z value Pr(>|z|)  
#> (Intercept)        -0.7954736  0.3088499  -2.576   0.0100 *
#> pb_ch               0.0143683  0.0681262   0.211   0.8330  
#> tb_ch               0.0917002  0.0496423   1.847   0.0647 .
#> bhs_ch              0.1202364  0.0771130   1.559   0.1189  
#> pb_ch:tb_ch        -0.0051317  0.0063136  -0.813   0.4163  
#> pb_ch:bhs_ch        0.0012289  0.0096951   0.127   0.8991  
#> tb_ch:bhs_ch       -0.0139993  0.0124736  -1.122   0.2617  
#> pb_ch:tb_ch:bhs_ch  0.0005919  0.0010741   0.551   0.5816  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
$\endgroup$

1 Answer 1

1
$\begingroup$

I want to use a mixed-effect model so that missing values can be handled by maximum-likelihood estimation (my assumption is that this is a feasible approach).

I think you may be misunderstanding the precise meaning of the assertion "mixed models can handle missing data well". "Modern" mixed models (such as lme4) can handle data where the groups are unbalanced, i.e. where different groups have different numbers of observations/observations corresponding to different subsets of fixed factors); this is in contrast to classic mixed-model ANOVA. Copying from my StackOverflow answer here:


In typical multivariate ANOVA contexts, the data are set up in wide format, i.e. one row for each group. For example:

id  t1  t2  t3  t4  t5
A   1.0 2.0 3.0 2.1 7.2
B   1.1 1.9 2.4 2.3 1.4
...

In this format, missing or unbalanced data would appear as an NA; if we didn't have an observation for group A at time t5, that row of the data would be 1.0 2.0 3.0 2.1 NA. However, in mixed-model-world we usually represent the data in long format:

id time value
A  t1   1.0
A  t2   2.0
A  t3   3.0
A  t4   2.1
...

so in this case we wouldn't even include the missing observation in the first place. In the MANOVA world we would have to decide how to deal with a group where any of the observations are missing; in the mixed-model world this corresponds to discarding a single group/time observation, corresponding to the complete-case analysis I described at the beginning.


In your case, you are missing values of responses and covariates, not observations within groups. (In your sample data, there is only one observation per group, which makes a mixed-model analysis inappropriate/unnecessary in any case ... maybe this doesn't represent a full data set? The typical solution to this case, if (1) your missing data are important enough to be worth the effort of including partial cases and (2) you can satisfy yourself that there are no important biases in the pattern of missingness (i.e. missing values are "missing completely at random" [MCAR] or "missing at random" [MAR], e.g. see here), is to do some kind of multiple imputation. In R, you can do this

Both will require some effort ...

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.