{"action":"created","issue":{"url":"https://api.github.com/repos/florianhartig/DHARMa/issues/101","repository_url":"https://api.github.com/repos/florianhartig/DHARMa","labels_url":"https://api.github.com/repos/florianhartig/DHARMa/issues/101/labels{/name}","comments_url":"https://api.github.com/repos/florianhartig/DHARMa/issues/101/comments","events_url":"https://api.github.com/repos/florianhartig/DHARMa/issues/101/events","html_url":"https://github.com/florianhartig/DHARMa/issues/101","id":392480573,"node_id":"MDU6SXNzdWUzOTI0ODA1NzM=","number":101,"title":"Missing At Random missing data","user":{"login":"drizopoulos","id":13236687,"node_id":"MDQ6VXNlcjEzMjM2Njg3","avatar_url":"https://avatars2.githubusercontent.com/u/13236687?v=4","gravatar_id":"","url":"https://api.github.com/users/drizopoulos","html_url":"https://github.com/drizopoulos","followers_url":"https://api.github.com/users/drizopoulos/followers","following_url":"https://api.github.com/users/drizopoulos/following{/other_user}","gists_url":"https://api.github.com/users/drizopoulos/gists{/gist_id}","starred_url":"https://api.github.com/users/drizopoulos/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/drizopoulos/subscriptions","organizations_url":"https://api.github.com/users/drizopoulos/orgs","repos_url":"https://api.github.com/users/drizopoulos/repos","events_url":"https://api.github.com/users/drizopoulos/events{/privacy}","received_events_url":"https://api.github.com/users/drizopoulos/received_events","type":"User","site_admin":false},"labels":[{"id":822784449,"node_id":"MDU6TGFiZWw4MjI3ODQ0NDk=","url":"https://api.github.com/repos/florianhartig/DHARMa/labels/Statistical%20question","name":"Statistical question","color":"5cd6d2","default":false},{"id":175711883,"node_id":"MDU6TGFiZWwxNzU3MTE4ODM=","url":"https://api.github.com/repos/florianhartig/DHARMa/labels/question","name":"question","color":"cc317c","default":true}],"state":"open","locked":false,"assignee":null,"assignees":[],"milestone":null,"comments":2,"created_at":"2018-12-19T07:57:02Z","updated_at":"2018-12-31T13:09:35Z","closed_at":null,"author_association":"NONE","body":"In the case of mixed models often we have to deal with incomplete data, for example, in longitudinal studies, subjects often drop out from the study. In this context, the framework of mixed models provides correct inferences under both the missing completely at random (missingness does not depend on the outcome) and missing at random (missingness depends on the observed outcomes) mechanisms.\r\n\r\nThe goodness-of-fit method implemented in **DHARMa** works in the missing completely at random case, but it does not seem to work under the missing at random one. A script illustrating the issue is provided below.\r\n\r\nMaking the connection to the posterior predictive tests, I think that a potential solution could be to impute the missing data from the model and apply procedure in the combined imputed and observed data.\r\n\r\n``` r\r\nlibrary(\"lme4\")\r\n#> Loading required package: Matrix\r\nlibrary(\"DHARMa\")\r\n#> Note that, since v0.1.6.2, DHARMa includes support for glmmTMB, but there are still a few minor limitations associatd with this package. Please see https://github.com/florianhartig/DHARMa/issues/16 for details, in particular if you use this for production.\r\n\r\nsimulate_data <- function (seed) {\r\n    set.seed(seed)\r\n    n <- 200 # number of subjects\r\n    K <- 8 # number of measurements per subject\r\n    t_max <- 5 # maximum follow-up time\r\n    \r\n    # we constuct a data frame with the design: \r\n    # everyone has a baseline measurment, and then measurements at random follow-up times\r\n    DF <- data.frame(id = rep(seq_len(n), each = K),\r\n                     time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),\r\n                     sex = rep(gl(2, n/2, labels = c(\"male\", \"female\")), each = K))\r\n    \r\n    # design matrices for the fixed and random effects\r\n    X <- model.matrix(~ sex * time, data = DF)\r\n    Z <- model.matrix(~ time, data = DF)\r\n    \r\n    betas <- c(12, -1.2, 1.6, -0.9) # fixed effects coefficients\r\n    sigma <- 0.5 # standard deviation error terms\r\n    D11 <- 0.5 # variance of random intercepts\r\n    D22 <- 0.1 # variance of random slopes\r\n    \r\n    # we simulate random effects\r\n    b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))\r\n    # linear predictor\r\n    eta_y <- as.vector(X %*% betas + rowSums(Z * b[DF$id, ]))\r\n    # we simulate normal longitudinal data\r\n    DF$y <- rnorm(n * K, mean = eta_y, sd = sigma)\r\n    \r\n    # missing completely at random missing data\r\n    DF$y2 <- DF$y # copy of y\r\n    set_to_na_mcar <- function (x) {\r\n        ind <- seq(2, length(x))\r\n        how_many_nas <- sample(K-2, 1)\r\n        x[sample(ind, how_many_nas)] <- NA\r\n        x\r\n    }\r\n    DF$y2 <- with(DF, ave(y2, id, FUN = set_to_na_mcar))\r\n    \r\n    # missing at random missing data\r\n    DF$y3 <- DF$y # copy of y\r\n    # if any y was above 13 except the baseline measurement, the subject drops out\r\n    set_to_na_mar <- function (x) {\r\n        nx <- length(x)\r\n        ind <- which(x[-1] > 13)[1] + 2\r\n        if (!is.na(ind) && ind < nx)\r\n            x[seq(ind, nx)] <- NA\r\n        x\r\n    }\r\n    DF$y3 <- with(DF, ave(y3, id, FUN = set_to_na_mar))\r\n    \r\n    DF\r\n}\r\n\r\n##########################################################################################\r\n##########################################################################################\r\n\r\nM <- 100\r\nbetas1 <- betas2 <- betas3 <- matrix(0.0, M, 4)\r\np_values1 <- p_values2 <- p_values3 <- numeric(M)\r\nfor (m in seq_len(M)) {\r\n    DF <- simulate_data(2018 + m)\r\n    ##\r\n    # fit the model and residuals in complete data\r\n    fm1 <- lmer(y ~ sex * time + (time | id), data = DF)\r\n    betas1[m, ] <- fixef(fm1)\r\n    res1 <- simulateResiduals(fittedModel = fm1, n = 500)\r\n    p_values1[m] <- testUniformity(simulationOutput = res1, plot = FALSE)$p.value\r\n    ###\r\n    # fit the model and residuals in missing completely at random\r\n    fm2 <- lmer(y2 ~ sex * time + (time | id), data = DF)\r\n    betas2[m, ] <- fixef(fm2)\r\n    res2 <- simulateResiduals(fittedModel = fm2, n = 500)\r\n    p_values2[m] <- testUniformity(simulationOutput = res2, plot = FALSE)$p.value\r\n    ###\r\n    # fit the model and residuals in missing at random\r\n    fm3 <- lmer(y3 ~ sex * time + (time | id), data = DF)\r\n    betas3[m, ] <- fixef(fm3)\r\n    res3 <- simulateResiduals(fittedModel = fm3, n = 500)\r\n    p_values3[m] <- testUniformity(simulationOutput = res3, plot = FALSE)$p.value\r\n}\r\n\r\n# Bias\r\ntrue_betas <- c(12, -1.2, 1.6, -0.9)\r\ncolMeans(betas1 - rep(true_betas, each = M))\r\n#> [1]  0.013800880 -0.010452142 -0.001517219 -0.002801698\r\ncolMeans(betas2 - rep(true_betas, each = M))\r\n#> [1]  0.0098002693 -0.0056694382  0.0002078754 -0.0050659936\r\ncolMeans(betas3 - rep(true_betas, each = M))\r\n#> [1]  0.011886233 -0.011053783  0.005843979 -0.007349007\r\n\r\n# p-value for Goodness-of-Fit\r\nsummary(p_values1)\r\n#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \r\n#> 0.06809 0.36906 0.63980 0.58937 0.78810 0.98741\r\nsummary(p_values2)\r\n#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \r\n#> 0.01241 0.43845 0.62246 0.60387 0.81227 0.99493\r\nsummary(p_values3)\r\n#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. \r\n#> 0.000e+00 3.800e-09 5.491e-08 6.299e-06 1.684e-06 1.888e-04\r\n```\r\n\r\n "},"comment":{"url":"https://api.github.com/repos/florianhartig/DHARMa/issues/comments/450644044","html_url":"https://github.com/florianhartig/DHARMa/issues/101#issuecomment-450644044","issue_url":"https://api.github.com/repos/florianhartig/DHARMa/issues/101","id":450644044,"node_id":"MDEyOklzc3VlQ29tbWVudDQ1MDY0NDA0NA==","user":{"login":"florianhartig","id":5457753,"node_id":"MDQ6VXNlcjU0NTc3NTM=","avatar_url":"https://avatars0.githubusercontent.com/u/5457753?v=4","gravatar_id":"","url":"https://api.github.com/users/florianhartig","html_url":"https://github.com/florianhartig","followers_url":"https://api.github.com/users/florianhartig/followers","following_url":"https://api.github.com/users/florianhartig/following{/other_user}","gists_url":"https://api.github.com/users/florianhartig/gists{/gist_id}","starred_url":"https://api.github.com/users/florianhartig/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/florianhartig/subscriptions","organizations_url":"https://api.github.com/users/florianhartig/orgs","repos_url":"https://api.github.com/users/florianhartig/repos","events_url":"https://api.github.com/users/florianhartig/events{/privacy}","received_events_url":"https://api.github.com/users/florianhartig/received_events","type":"User","site_admin":false},"created_at":"2018-12-31T13:09:35Z","updated_at":"2018-12-31T13:09:35Z","author_association":"OWNER","body":"Hi Dimitris,\r\n\r\nI guess I could implement something like this, but the question is how DHARMa should know when to apply this procedure. Moreover, there is the MAR case, but there are probably many similar cases in other fields. \r\n\r\nBut even if we take only the MAR: how should DHARMa actually know that data is MAR? From what I understand, it's not trivial to check this. It seems to me that if a user has sufficient statistical proficiency to be sure that data is MAR, they should also be proficient enough to understand how to extract residuals in DHARMa and do appropriate corrections. \r\n\r\nWhat I could possibly do is to create a section in the vignette about \"apparent\" residual problems, and mentioning the MAR as one case (I can think about a few other situations where residuals are easily misinterpreted, both regarding false positives and false negatives)\r\n\r\n"}}
{"actor":{"display_login":"florianhartig"}}