Bayesian Ordered Logit - Versuch, vorhergesagte y im Zeitverlauf basierend auf einem modalen Individuum zu berechnen

Question

Ich habe einen Datensatz, der mehrere Umfragen mit verschiedenen Ländern über Jahre hinweg kombiniert. Meine abhängige Variable (lrparty) ist die ideologische Position einer Partei (von 0 bis 10) nach Umfrageteilnehmern. Ich habe mehrere unabhängige Variablen wie Alter, Geschlecht, Bildung, Parteinahme und Einkommen der Befragten.

Dann würde ich für jede Partei und jede Umfrage die vorhergesagten Werte von lrparty gemäß dem modalen Individuum planen (zB Befragter mit Alter = 31, weiblich = 1, Bildung = 2, Einkommen = 2, und Partisan = 1) im Laufe der Zeit. Der Graph würde also wie folgt aussehen: x-Achse = Jahre; y-Achse = vorhergesagte Werte von lrparty gemäß dem modalen Individuum.

In Summe sind dies die Stufen von dem, was ich zu tun versuchen: 1. Schätzung des Modells: bestellt logit der Platzierung der Partei (lrparty), des Geschlechts regredieren, Alter, Bildung, Einkommen und Parteilichkeit von Befragten.

2. Hintere Züge nehmen.

3. Predict Platzierung der Partei für die modale Befragten (zB 500 Unentschieden)

4. Dann erwarte ich einen Datensatz zu haben, die wie folgt aussehen sollte: Jahr, Umfrage, Land, Party (cmp Code),% fehlende Platzierungen, x1: x500 (aus den Zeichnungen)

5. Aus diesem Datensatz würde ich meine Plots generieren. Zum Beispiel für UK, nach Umfrage CSES.

Um es herauszufinden, den Code aus, begann ich nur eine Umfrage (CSE), ein Land (uk), und eine Partei (Konservative) mit, wie Sie unten in meinem Code zu sehen. Aber ich weiß nicht, wie ich von der Stelle, an der ich bin, in den Code zu der Handlung komme, die ich will (oben beschrieben).

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) 
    ##? 

Vielen Dank im Voraus!

Answer 1

Ok. Nach mehreren Versuchen versuchte ich den Code. Da es für jemand anderen von Interesse sein könnte, poste ich den Code unten.

## 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))