Bayesian Ordered Logit - Essayer de tracer prédite y au fil du temps basé sur un individu modal

Question

J'ai un jeu de données qui combine plusieurs enquêtes, avec différents pays, au fil des années. Ma variable dépendante (lrparty) est la position idéologique d'une partie (allant de 0 à 10) selon les répondants à l'enquête. J'ai plusieurs variables indépendantes telles que l'âge, le sexe, l'éducation, la partisanerie et le revenu des répondants. Puis, pour chaque parti et chaque enquête, je voudrais tracer les valeurs prédites de lrparty selon l'individu modal (par exemple, répondant avec l'âge = 31, femelle = 1, l'éducation = 2, le revenu = 2, et partisan = 1) au fil du temps. Donc, le graphique ressemblerait à: x-axis = years; y-axis = valeurs prédites de lrparty en fonction de l'individu modal.

En somme, ce sont les étapes de ce que je suis en train de faire: 1. Estimation du modèle: logit ordonné de placement de la partie (lrparty), régressant sur le sexe, l'âge, l'éducation, le revenu et la partisanerie des répondants.

2. Prendre des tirages postérieurs.

3. Prédire le placement de la partie pour l'intimée modal (par exemple, 500 tirages)

4. Puis, je me attends d'avoir un ensemble de données qui devrait ressembler à: Année, Enquête, pays, Parti (code cmp),% emplacements manquants, x1: x500 (à partir des tirages)

5. À partir de cet ensemble de données, je générerais mes graphiques. Par exemple, pour le Royaume-Uni, selon l'enquête CSES.

Afin de comprendre le code, je commencé à utiliser une seule enquête (CSEs), un pays (Royaume-Uni), et une partie (conservateurs) que vous pouvez voir dans mon code ci-dessous. Mais je ne sais pas comment aller d'où je suis dans le code à l'intrigue que je veux (ci-dessus décrit).

library(rstan) 
    library(tidyverse) 
    library(brms) 
    library(ggplot2) 
    library(ggthemes) 
    library(ggmcmc) 

    ## Data: 
    load("pbrands.RData") 

    ## Keeping only country = uk; survey = cses; party = conservatives 
    uk_cses_con = pbrands %>% 
    select(lrparty, female, age, education, income, partisan, year, survey,         
    country, cmp, party_name_short, party_name_english, lrs) %>% 
    filter(survey == "cses") %>% 
    filter(country == "uk") %>% 
    filter(cmp == 51620) 

    ## Conducting a Bayesian ordered logit model 
    fit <- brm(lrparty ~ age + income + education + female + partisan, 
     data = uk_cses_con, family = "cumulative", chains = 4, iter = 1000) 

    ## Trace and Density Plots for MCMC Samples 
    plot(fit) 

    ## Posterior Predictive Checks 
    pp_check(fit) 

    ## Getting variables' modes: 
    getmode <- function(v) { 
    uniqv <- unique(v) 
    uniqv[which.max(tabulate(match(v, uniqv)))] 
    } 

    getmode(uk_cses_con$age) 
    getmode(uk_cses_con$female) 
    getmode(uk_cses_con$education) 
    getmode(uk_cses_con$income) 
    getmode(uk_cses_con$partisan) 

    ## Creating the data frame for the modal individual 
    newavg <- data.frame(age = 31, female = 1, education = 2, income = 2,    
    partisan = 0, years = uk_cses_con$year) 

    ## predict response for new data 
    pred <- predict(fit, newdata = newavg) 

    # extract posterior samples of population-level effects 
    samples1 <- posterior_samples(fit) 

    ## Display marginal effects of predictors 
    marginal <- marginal_effects(fit) 

    ## Plot predicted lrparty (my dependent variable) over time (with error:   
    confidence interval) based on the modal respondent (age = 31, female = 1,     
    education = 0, income = 0, partisan = 0) 
    ##? 

Merci d'avance!

Answer 1

Ok. Après plusieurs tentatives d'essais et d'erreurs, j'ai compris le code. Comme cela pourrait intéresser quelqu'un d'autre, je poste le code ci-dessous.

## Packages 
    install.packages(c("bmrs", "coda", "mvtnorm", "devtools")) 
    library(devtools) 
    library(tidyverse) 
    library(brms) 

    ## Loading the data 
    load('~/Data/mydata.RData') 

    ## Keeping the variables of our interest 
    mydata = mydata %>% 
    select(lrparty, female, age, education, income, partisan, year, survey, 
    country, cmp, party_name_short, party_name_english, lrs) 

    ## Function for calculating modes 
    getmode <- function(v) { 
    uniqv = unique(v) 
    uniqv[which.max(tabulate(match(v, uniqv)))] 
    } 

    ## Finding Modal respondents by country, survey, and party: 

    ## Modes by country 
    mode_by_country = mydata %>% 
    group_by(country) %>% 
    mutate(modal_age = getmode(na.omit(age))) %>% 
    mutate(modal_inc = getmode(na.omit(income))) %>% 
    mutate(modal_female = getmode(na.omit(female))) %>% 
    mutate(modal_edu = getmode(na.omit(education))) %>% 
    mutate(modal_partisan = getmode(na.omit(partisan))) %>% 
    filter(!duplicated(country)) 

    mode_by_country = mode_by_country[ , c(9, 14:18)] 

    mode_by_country = mode_by_country[-40, ] 

    ## Build object to receive the information we want to store 
    runner <- c() 
    pred = matrix(NA, 2000, 11) 
    yhat = matrix(NA, 2000, 1) 

    ###### Conducting the model for UK with two parties ######### 
    uk = mydata %>% 
    select(lrparty, female, age, education, income, partisan, year, survey,    
    country, cmp, party_name_short, party_name_english, lrs) %>% 
    filter(survey == "cses") %>% 
    filter(country == "uk") %>% 
    filter(cmp == 51320 | cmp == 51620) 

    ## Finding how many regressions will be conducted 
    reglength <- length(unique(uk$survey)) * length(unique(uk$year)) * length(unique(uk$cmp)) 

    ## Creating our modal British individual based on mode_by_country 
    mode_by_country[mode_by_country$country == "uk", c(2:6)] 

    newavg <- data.frame(age = 35, income = 2, female = 1, education = 2, partisan = 0) 

    ## Loop to conduct the ordered logit in rstan, using iter=1000, and   chains=4 

    for(p in na.omit(unique(uk$cmp))){ 

    hold <- uk[uk$cmp == p, ] 
    country <- hold$country[1] 

    for(s in na.omit(unique(hold$survey))){ 

     hold1 <- hold[hold$survey == s, ] 

     for(y in na.omit(unique(hold1$year))){ 

      mod <- brm(lrparty ~ age + female + education + income + partisan, data = hold1[hold1$year == y, ], family = "cumulative", chains = 4, iter = 1000) 

      for(i in 1:2000) { 
       pred[i,] <- predict(mod, newdata = newavg, probs = c(0.025, 0.975), summary=TRUE) 
       yhat[i] <- sum(pred[i, ] * c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)) 
      } 

       newData <- data.frame(country, p, s, y, pred, yhat) 
       newData$m <- mean(newData$yhat) 
       newData$sd <- sd(newData$yhat) 
       newData$lower <- newData$m - 1.96*newData$sd 
       newData$upper <- newData$m + 1.96*newData$sd 

      runner <- rbind(runner, newData) 
     } 
     } 
    } 

    ## Keeping unique values within dataset 
    uniqdata = runner %>% 
    filter(!duplicated(m)) 

    ## Creating the Figure 
    uniqdata2 <- uniqdata[, c(1:4, 17:20)] 

    uniqdata3 <- uniqdata2 %>% 
    gather(variable, value, -(y:p)) %>% 
    unite(temp, p, variable) %>% 
    spread(temp, value) 

    uniqdata3 = uniqdata3[ , -c(3,6,8,11)] 

    names(uniqdata3)[3:8] = c("lower_lab", "m_lab", "upper_lab", "lower_con", "m_con", "upper_con") 

    uniqdata3[3:8] = as.numeric(unlist(uniqdata3[3:8])) 

    ## Plot: Predicted Party Ideological Placement for Modal British Respondent 

    ggplot(uniqdata3, aes(x = (y))) + geom_line(aes(y = m_lab, colour = "Labor")) + geom_ribbon(aes(ymin = lower_lab,ymax = upper_lab, 
       linetype=NA), alpha = .25) + 
    geom_line(aes(y = m_con, color = "Conservatives")) + 
    geom_ribbon(aes(ymin = lower_con, 
       ymax = upper_con, 
       linetype=NA), alpha = .25) + 
    theme_bw() + 
    theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5)) + labs(title = "Predicted Party Ideological Placement for Modal British Respondent \n Survey = CSES") + theme(plot.title = element_text(color="black", size=20, face="bold.italic"), axis.title.x = element_text(color="black", size=15, face="italic"), axis.title.y = element_text(color="black", size=15, face="italic")) + 
    theme(legend.title = element_blank()) + 
    theme(axis.text.x = element_text(color="black", size= 12.5), axis.text.y = element_text(color="black", size=12.5)) + theme(legend.text = element_text(size=15)) + scale_x_continuous(name="Year", breaks=seq(1997, 2005, 2)) + scale_y_continuous(name="Left-Right Party Position", limits=c(0, 10))