Abstract

Expanding upon Blosnich, Shipherd, and Kauth’s seminal 2019 research, this study investigates anti-gay beliefs among veterans versus non-veterans using updated 2018 data from the General Social Survey (GSS). Analyzing GSS data from 2010 to 2016, the original study found a notable correlation between veteran status and increased disapproval of same-sex marriage. However, it also noted a declining trend in anti-gay sentiments among veterans across successive surveys, indicating a gradual convergence with the more tolerant attitudes of non-veterans over time. Building on this, the current study delves further into these trends by examining post-2015 Supreme Court decision data, which legalized same-sex marriage nationwide—a decision emblematic of the civil rights upheld by the military. By integrating these additional GSS waves, the study aims to offer deeper insights into veterans’ evolving attitudes about broader societal acceptance shifts towards same-sex marriage. Results of the chi-square tests of independence reveal significant associations between veteran status and demographic variables such as sex and race. Specifically, males and whites are more likely to be veterans, highlighting potential demographic disparities within military service. Furthermore, attitudes towards speaking about, teaching, and the legality of homosexuality (SPKHOMO, COLHOMO, and LIBHOMO) do not show a significant difference in opinion between veterans and non-veterans, suggesting that military service does not markedly influence these views. However, the regression analysis of attitudes towards homosexuality (HOMOSEX) and same-sex marriage (MARHOMO), when quantified numerically, unveils nuanced differences. Veterans, on average, score lower on the HOMOSEX scale, indicating a tendency towards less acceptance. Similarly, veterans expressed lower towards same-sex marriage (MARHOMO). Interestingly, while non-veterans consistently harbored anti-gay sentiments towards same-sex marriage, trends from 2010 to 2018 show a shift among veterans towards more positive or accepting viewpoints of same-sex marriage. In essence, this study sheds light on the evolving attitudes of veterans about societal acceptance shifts towards same-sex marriage and anti-gay beliefs, highlighting a potential trend towards greater tolerance among this demographic group in recent years.

# Install and load required packages
if (!requireNamespace("gssr", quietly = TRUE)) {
  remotes::install_github("kjhealy/gssr", dependencies = TRUE)
}
library(gssr)
## Package loaded. To attach the GSS data, type data(gss_all) at the console.
## For the codebook, type data(gss_dict).
## For the panel data and documentation, type e.g. data(gss_panel08_long) and data(gss_panel_doc).
## For help on a specific GSS variable, type ?varname at the console.
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(explore)
library(haven)
gss_all <- read_sav("C:/Users/adolp/OneDrive - Significant Results/PostDoc/Lynne/Military/GSS7218_R3.sav")
View(gss_all)
library(survey)
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
library(dplyr)

# Ensure 'YEAR' is numeric
gss_all$YEAR <- as.numeric(as.character(gss_all$YEAR))

# Filter the dataset for specific years
gss_all_filtered <- gss_all %>%
  filter(YEAR %in% c(2010, 2012, 2014, 2016, 2018, 2020, 2022)) %>%
  # Selectively retain only the listed variables
  select(YEAR, VETYEARS, SPKHOMO, COLHOMO, LIBHOMO, HOMOSEX, MARHOMO, 
         SEX, RACE, AGE, EDUC, SEXORNT, RELITEN, WTSSNR, VSTRAT, VPSU)
# Recode variables

# Veteran status
library(dplyr)

# Assuming gss_all_filtered is already loaded and exists
gss_all_filtered <- gss_all_filtered %>%
  mutate(Veteran_Status = ifelse(VETYEARS %in% c(2, 3, 4), "Veteran", "Non-Veteran"))

# Anti-gay perceptions
levels_speaking <- c("Allowed to speak", "Not allowed")
levels_teaching <- c("Allowed to teach", "Not allowed")
levels_removing <- c("Remove", "Not remove")
levels_sex_relations <- c("Always wrong", "Almost always wrong", "Sometimes wrong", "Not wrong at all")
levels_marriage <- c("Strongly agree", "Agree", "Neither agree nor disagree", "Disagree", "Strongly disagree")

# Recode and label the variables
gss_all_filtered <- gss_all_filtered %>%
  mutate(
    SPKHOMO = factor(SPKHOMO, levels = c(1, 2), labels = levels_speaking),
    COLHOMO = factor(COLHOMO, levels = c(4, 5), labels = levels_teaching),
    LIBHOMO = factor(LIBHOMO, levels = c(1, 2), labels = levels_removing),
    HOMOSEX = factor(HOMOSEX, levels = c(1, 2, 3, 4), labels = levels_sex_relations),
    MARHOMO = factor(MARHOMO, levels = c(1, 2, 3, 4, 5), labels = levels_marriage)
  )


# Demographics
gss_all_filtered <- gss_all_filtered %>%
  mutate(
    SEX = factor(SEX, levels = c(1, 2), labels = c("Male", "Female")),
    RACE = factor(RACE, levels = c(1, 2, 3), labels = c("White", "Black", "Other")),
    SEXORNT = factor(SEXORNT, levels = c(1, 2, 3), labels = c("Gay, Lesbian, or Homosexual", "Bisexual", "Heterosexual or Straight")),
    RELITEN = factor(RELITEN, levels = c(1, 2, 3, 4), labels = c("Strong", "Not very strong", "Somewhat strong (vol.)", "No religion"))
  )
# Load the survey package if not already loaded
library(survey, quietly = T)

# Create survey design object
gss_all_filtered <- svydesign(
  ids = ~VPSU,
  strata = ~VSTRAT,
  weights = ~WTSSNR,
  data = gss_all_filtered,
  nest = TRUE
)

# View the survey design object
gss_all_filtered
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (816) clusters.
## svydesign(ids = ~VPSU, strata = ~VSTRAT, weights = ~WTSSNR, data = gss_all_filtered, 
##     nest = TRUE)
# Listing variables using names on the data component
names(gss_all_filtered$variables)
##  [1] "YEAR"           "VETYEARS"       "SPKHOMO"        "COLHOMO"       
##  [5] "LIBHOMO"        "HOMOSEX"        "MARHOMO"        "SEX"           
##  [9] "RACE"           "AGE"            "EDUC"           "SEXORNT"       
## [13] "RELITEN"        "WTSSNR"         "VSTRAT"         "VPSU"          
## [17] "Veteran_Status"
# Alternatively, you can directly access the data frame
names(gss_all_filtered$data)
## NULL
# Function to calculate and print the demographic breakdown by Veteran Status
print_demographic_breakdown <- function(variable, design) {
  # Construct the formula strings
  table_formula_str <- paste0("~", variable, "+ Veteran_Status")
  chi_formula_str <- paste0("~", variable, "+ Veteran_Status")
  
  # Convert string to formula
  table_formula <- as.formula(table_formula_str)
  chi_formula <- as.formula(chi_formula_str)
  
  # Calculate and print the demographic table
  demographic_table <- svytable(table_formula, design = design)
  cat("Demographic table for", variable, "by Veteran Status:\n")
  print(demographic_table)
  
  # Perform and print the Chi-Square test of independence
  chi_test_result <- svychisq(chi_formula, design = design)
  cat("Chi-Square Test of Independence for", variable, "and Veteran Status:\n")
  print(chi_test_result)
  
  # If you want to calculate and print means for each category (for ordinal variables), you can adapt this part
  # Note: Direct mean calculation on factors will not be meaningful without converting them to numeric if they represent ordinal scales
}

# Assuming 'gss_all_filtered' is your survey design object
# Example usage for SEX demographic
print_demographic_breakdown("SEX", gss_all_filtered)
## Demographic table for SEX by Veteran Status:
##         Veteran_Status
## SEX      Non-Veteran    Veteran
##   Male    4550.60843  817.31746
##   Female  6309.56923   94.50489
## Chi-Square Test of Independence for SEX and Veteran Status:
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(chi_formula, design = design)
## F = 651.3, ndf = 1, ddf = 408, p-value < 2.2e-16
# Analyze differences in RACE by Veteran Status
print_demographic_breakdown("RACE", gss_all_filtered)
## Demographic table for RACE by Veteran Status:
##        Veteran_Status
## RACE    Non-Veteran    Veteran
##   White  7929.33147  721.74014
##   Black  1649.28119  127.81512
##   Other  1281.56500   62.26708
## Chi-Square Test of Independence for RACE and Veteran Status:
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(chi_formula, design = design)
## F = 8.4844, ndf = 1.9867, ddf = 810.5534, p-value = 0.000234
# Analyze differences in SEXORNT by Veteran Status
print_demographic_breakdown("SEXORNT", gss_all_filtered)
## Demographic table for SEXORNT by Veteran Status:
##                              Veteran_Status
## SEXORNT                       Non-Veteran     Veteran
##   Gay, Lesbian, or Homosexual  147.833157    8.783141
##   Bisexual                     206.836954   11.599408
##   Heterosexual or Straight    7929.019039  680.727835
## Chi-Square Test of Independence for SEXORNT and Veteran Status:
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(chi_formula, design = design)
## F = 1.4498, ndf = 1.8845, ddf = 768.8797, p-value = 0.2357
# Analyze differences in RELITEN by Veteran Status
print_demographic_breakdown("RELITEN", gss_all_filtered)
## Demographic table for RELITEN by Veteran Status:
##                         Veteran_Status
## RELITEN                  Non-Veteran    Veteran
##   Strong                  3827.88887  316.42818
##   Not very strong         3865.56290  364.10489
##   Somewhat strong (vol.)   688.93587   47.59888
##   No religion             2279.24936  166.40173
## Chi-Square Test of Independence for RELITEN and Veteran Status:
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(chi_formula, design = design)
## F = 2.673, ndf = 2.9676, ddf = 1210.7659, p-value = 0.04671
# Load necessary packages for data manipulation and survey analysis
library(dplyr)
library(survey)

# Assuming 'gss_all_filtered' is your survey design object created using svydesign

# Crosstab and Chi-square test for Veteran_Status and SPKHOMO using survey package
svy_crosstab_spkhomo <- svytable(~Veteran_Status + SPKHOMO, design = gss_all_filtered)
svy_row_pct_spkhomo <- prop.table(svy_crosstab_spkhomo, margin = 1) * 100

# Chi-square test is not directly applicable in the same way with survey objects,
# but you can use svychisq() for a design-based chi-square test as follows:
svy_chisq_test_spkhomo <- svychisq(~Veteran_Status + SPKHOMO, design = gss_all_filtered)

# Print the crosstab with row percentages
cat("Crosstab for Veteran_Status and SPKHOMO:\n")
## Crosstab for Veteran_Status and SPKHOMO:
print(svy_crosstab_spkhomo)
##               SPKHOMO
## Veteran_Status Allowed to speak Not allowed
##    Non-Veteran       6261.70827   787.53566
##    Veteran            545.25592    71.63946
cat("Row Percentages for Veteran_Status and SPKHOMO:\n")
## Row Percentages for Veteran_Status and SPKHOMO:
print(svy_row_pct_spkhomo)
##               SPKHOMO
## Veteran_Status Allowed to speak Not allowed
##    Non-Veteran         88.82808    11.17192
##    Veteran             88.38710    11.61290
cat("Chi-square Test for Veteran_Status and SPKHOMO:\n")
## Chi-square Test for Veteran_Status and SPKHOMO:
print(svy_chisq_test_spkhomo)
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~Veteran_Status + SPKHOMO, design = gss_all_filtered)
## F = 0.096322, ndf = 1, ddf = 408, p-value = 0.7564
# Repeat the process for Veteran_Status and COLHOMO
svy_crosstab_colhomo <- svytable(~Veteran_Status + COLHOMO, design = gss_all_filtered)
svy_row_pct_colhomo <- prop.table(svy_crosstab_colhomo, margin = 1) * 100
svy_chisq_test_colhomo <- svychisq(~Veteran_Status + COLHOMO, design = gss_all_filtered)

cat("Crosstab for Veteran_Status and COLHOMO:\n")
## Crosstab for Veteran_Status and COLHOMO:
print(svy_crosstab_colhomo)
##               COLHOMO
## Veteran_Status Allowed to teach Not allowed
##    Non-Veteran        6185.6986    841.1106
##    Veteran             509.6340    106.3300
cat("Row Percentages for Veteran_Status and COLHOMO:\n")
## Row Percentages for Veteran_Status and COLHOMO:
print(svy_row_pct_colhomo)
##               COLHOMO
## Veteran_Status Allowed to teach Not allowed
##    Non-Veteran         88.02998    11.97002
##    Veteran             82.73763    17.26237
cat("Chi-square Test for Veteran_Status and COLHOMO:\n")
## Chi-square Test for Veteran_Status and COLHOMO:
print(svy_chisq_test_colhomo)
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~Veteran_Status + COLHOMO, design = gss_all_filtered)
## F = 12.556, ndf = 1, ddf = 408, p-value = 0.0004406
# And for Veteran_Status and LIBHOMO
svy_crosstab_libhomo <- svytable(~Veteran_Status + LIBHOMO, design = gss_all_filtered)
svy_row_pct_libhomo <- prop.table(svy_crosstab_libhomo, margin = 1) * 100
svy_chisq_test_libhomo <- svychisq(~Veteran_Status + LIBHOMO, design = gss_all_filtered)

cat("Crosstab for Veteran_Status and LIBHOMO:\n")
## Crosstab for Veteran_Status and LIBHOMO:
print(svy_crosstab_libhomo)
##               LIBHOMO
## Veteran_Status    Remove Not remove
##    Non-Veteran 1229.7314  5800.9766
##    Veteran      123.1137   494.3051
cat("Row Percentages for Veteran_Status and LIBHOMO:\n")
## Row Percentages for Veteran_Status and LIBHOMO:
print(svy_row_pct_libhomo)
##               LIBHOMO
## Veteran_Status   Remove Not remove
##    Non-Veteran 17.49086   82.50914
##    Veteran     19.94006   80.05994
cat("Chi-square Test for Veteran_Status and LIBHOMO:\n")
## Chi-square Test for Veteran_Status and LIBHOMO:
print(svy_chisq_test_libhomo)
## 
##  Pearson's X^2: Rao & Scott adjustment
## 
## data:  svychisq(~Veteran_Status + LIBHOMO, design = gss_all_filtered)
## F = 1.9807, ndf = 1, ddf = 408, p-value = 0.1601
# Load necessary packages for data manipulation and survey analysis
library(dplyr)
library(survey)

# Assuming 'gss_all_filtered' is your survey design object created using svydesign

# Step 1: Ensure HOMOSEX is an ordered factor (adjust the levels according to your actual data)
gss_all_filtered$data$HOMOSEX <- factor(gss_all_filtered$data$HOMOSEX,
                                        levels = c("Always wrong", "Almost always wrong", "Sometimes wrong", "Not wrong at all"),
                                        ordered = TRUE)

# Step 2: Convert the ordered factor to numeric
# Since survey design objects don't directly support transformation, perform this on the data frame first
gss_all_filtered$data$NUMERIC_HOMOSEX <- as.numeric(gss_all_filtered$data$HOMOSEX)

# Update the survey design object to include the new numeric variable
gss_all_filtered <- update(gss_all_filtered, NUMERIC_HOMOSEX = as.numeric(HOMOSEX))

# Now, fit a linear model using svyglm for comparing NUMERIC_HOMOSEX scores between Veteran and Non-Veteran groups
model_homosex <- svyglm(NUMERIC_HOMOSEX ~ Veteran_Status, design = gss_all_filtered)

# To interpret the model in a manner similar to a t-test, you would look at the coefficient for Veteran_Status
# This gives you the difference in means between the two groups, accounting for the survey design

# Display the summary of the model to see the coefficients and p-values
summary(model_homosex)
## 
## Call:
## svyglm(formula = NUMERIC_HOMOSEX ~ Veteran_Status, design = gss_all_filtered)
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_HOMOSEX = as.numeric(HOMOSEX))
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.69348    0.02164 124.464  < 2e-16 ***
## Veteran_StatusVeteran -0.34951    0.06454  -5.416 1.05e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.010888)
## 
## Number of Fisher Scoring iterations: 2
# Load necessary packages for data manipulation and survey analysis
library(dplyr)
library(survey)

# Assuming 'gss_all_filtered' is your survey design object created using svydesign

# Step 1: Ensure MARHOMO is an ordered factor (adjust levels according to your actual data)
gss_all_filtered$data$MARHOMO <- factor(gss_all_filtered$data$MARHOMO,
                                        levels = c("Strongly agree", "Agree", "Neither agree nor disagree", "Disagree", "Strongly disagree"),
                                        ordered = TRUE)

# Step 2: Convert the ordered factor MARHOMO to numeric
# Perform this on the data frame first, as direct modifications are not straightforward with survey design objects
gss_all_filtered$data$NUMERIC_MARHOMO <- as.numeric(gss_all_filtered$data$MARHOMO)

# Update the survey design object to include the new numeric variable
gss_all_filtered <- update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))

# Now, fit a linear model using svyglm for comparing NUMERIC_MARHOMO scores between Veteran and Non-Veteran groups
model_marhomo <- svyglm(NUMERIC_MARHOMO ~ Veteran_Status, design = gss_all_filtered)

# To interpret the model in a manner similar to a t-test, look at the coefficient for Veteran_Status
# This provides the difference in means between the two groups, taking into account the survey design

# Display the summary of the model to see the coefficients and p-values
summary(model_marhomo)
## 
## Call:
## svyglm(formula = NUMERIC_MARHOMO ~ Veteran_Status, design = gss_all_filtered)
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            2.57936    0.02188 117.900  < 2e-16 ***
## Veteran_StatusVeteran  0.47631    0.06855   6.949  1.5e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 2.25076)
## 
## Number of Fisher Scoring iterations: 2
# Load necessary packages for data manipulation, survey analysis, and plotting
library(ggplot2)
library(dplyr)
library(survey)

# Assuming 'gss_all_filtered' is your survey design object created using svydesign
# and NUMERIC_MARHOMO is defined within this survey design object

# Calculate average MARHOMO by YEAR and Veteran_Status using survey package functions
avg_marhomo_by_year <- svyby(~NUMERIC_MARHOMO, by = ~YEAR + Veteran_Status, design = gss_all_filtered,
                             FUN = svymean, na.rm = TRUE)

# Transform the resulting object into a data frame for ggplot
avg_marhomo_by_year_df <- as.data.frame(avg_marhomo_by_year)
# Rename columns appropriately if needed
names(avg_marhomo_by_year_df)[names(avg_marhomo_by_year_df) == "NUMERIC_MARHOMO"] <- "Avg_MARHOMO"

# Plot
ggplot(avg_marhomo_by_year_df, aes(x = YEAR, y = Avg_MARHOMO, group = Veteran_Status, color = Veteran_Status)) +
  geom_line() + # For line plot
  geom_point() + # To add points to the line plot
  labs(title = "Average Level of Agreement with Same-Sex Marriage Across Years",
       x = "Year",
       y = "Average Agreement Level",
       caption = "Scale: 1=Strongly Agree to 5=Strongly Disagree") +
  theme_minimal() +
  scale_color_manual(values = c("Veteran" = "blue", "Non-Veteran" = "red")) # Customize colors

library(survey)

# Ensure categorical variables are factors with appropriate levels
# Example for RACE and SEXORNT, repeat as necessary for other categorical variables
gss_all_filtered$data$RACE <- factor(gss_all_filtered$data$RACE, levels = c(1, 2, 3), labels = c("White", "Black", "Other"))
gss_all_filtered$data$SEXORNT <- factor(gss_all_filtered$data$SEXORNT, levels = c(1, 2, 3), labels = c("Gay, Lesbian, or Homosexual", "Bisexual", "Heterosexual or Straight"))
gss_all_filtered$data$Veteran_Status <- factor(gss_all_filtered$data$Veteran_Status, levels = c("Non-Veteran", "Veteran"))

# SPKHOMO should also be a factor if it's categorical
gss_all_filtered$data$SPKHOMO <- factor(gss_all_filtered$data$SPKHOMO)

# Define the regression formula
formula <- as.formula("SPKHOMO ~ YEAR + RACE + AGE + EDUC + SEXORNT + Veteran_Status")

# Run the regression using svyglm
svyglm_result <- svyglm(formula, design = gss_all_filtered, family = quasibinomial())

# Print the summary of the regression result
summary(svyglm_result)
## 
## Call:
## svyglm(formula = formula, design = gss_all_filtered, family = quasibinomial())
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     117.271114  42.576554   2.754 0.006163 ** 
## YEAR                             -0.058586   0.021169  -2.768 0.005924 ** 
## RACEBlack                         0.768770   0.128713   5.973 5.35e-09 ***
## RACEOther                         0.470752   0.242926   1.938 0.053381 .  
## AGE                               0.015180   0.003903   3.890 0.000118 ***
## EDUC                             -0.217684   0.018278 -11.910  < 2e-16 ***
## SEXORNTBisexual                  -0.579794   0.759161  -0.764 0.445501    
## SEXORNTHeterosexual or Straight   0.385256   0.505333   0.762 0.446304    
## Veteran_StatusVeteran            -0.088361   0.192798  -0.458 0.646990    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.9573913)
## 
## Number of Fisher Scoring iterations: 6
library(survey)

# Ensure categorical variables are factors with appropriate levels
gss_all_filtered$data$RACE <- factor(gss_all_filtered$data$RACE, levels = c(1, 2, 3), labels = c("White", "Black", "Other"))
gss_all_filtered$data$SEXORNT <- factor(gss_all_filtered$data$SEXORNT, levels = c(1, 2, 3), labels = c("Gay, Lesbian, or Homosexual", "Bisexual", "Heterosexual or Straight"))
gss_all_filtered$data$Veteran_Status <- factor(gss_all_filtered$data$Veteran_Status)

# COLHOMO should also be a factor if it's categorical
gss_all_filtered$data$COLHOMO <- factor(gss_all_filtered$data$COLHOMO)

# Define the regression formula with COLHOMO as the DV
formula_colhomo <- as.formula("COLHOMO ~ YEAR + RACE + AGE + EDUC + SEXORNT + Veteran_Status")

# Run the regression using svyglm
# Adjust the family parameter based on the distribution of COLHOMO
# Assuming COLHOMO is binary or categorical with more than two levels, you might use binomial or multinomial
svyglm_result_colhomo <- svyglm(formula_colhomo, design = gss_all_filtered, family = quasibinomial())

# Print the summary of the regression result for COLHOMO
summary(svyglm_result_colhomo)
## 
## Call:
## svyglm(formula = formula_colhomo, design = gss_all_filtered, 
##     family = quasibinomial())
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     146.992316  43.474432   3.381 0.000797 ***
## YEAR                             -0.073494   0.021573  -3.407 0.000727 ***
## RACEBlack                         0.622002   0.122136   5.093 5.56e-07 ***
## RACEOther                         0.342151   0.242931   1.408 0.159820    
## AGE                               0.023690   0.003534   6.703 7.39e-11 ***
## EDUC                             -0.211966   0.017294 -12.257  < 2e-16 ***
## SEXORNTBisexual                  -0.388045   0.658566  -0.589 0.556057    
## SEXORNTHeterosexual or Straight   0.351869   0.463548   0.759 0.448274    
## Veteran_StatusVeteran             0.288797   0.167296   1.726 0.085111 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 0.9344709)
## 
## Number of Fisher Scoring iterations: 5
library(survey)

# Ensure categorical variables are factors with appropriate levels
gss_all_filtered$data$RACE <- factor(gss_all_filtered$data$RACE, levels = c(1, 2, 3), labels = c("White", "Black", "Other"))
gss_all_filtered$data$SEXORNT <- factor(gss_all_filtered$data$SEXORNT, levels = c(1, 2, 3), labels = c("Gay, Lesbian, or Homosexual", "Bisexual", "Heterosexual or Straight"))
gss_all_filtered$data$Veteran_Status <- factor(gss_all_filtered$data$Veteran_Status)

# Ensure LIBHOMO is a factor if it's categorical
gss_all_filtered$data$LIBHOMO <- factor(gss_all_filtered$data$LIBHOMO)

# Define the regression formula with LIBHOMO as the DV
formula_libhomo <- as.formula("LIBHOMO ~ YEAR + RACE + AGE + EDUC + SEXORNT + Veteran_Status")

# Run the regression using svyglm
# Adjust the family parameter based on the distribution of LIBHOMO
# Assuming LIBHOMO is binary or categorical with more than two levels, you might use binomial or multinomial
svyglm_result_libhomo <- svyglm(formula_libhomo, design = gss_all_filtered, family = quasibinomial())

# Print the summary of the regression result for LIBHOMO
summary(svyglm_result_libhomo)
## 
## Call:
## svyglm(formula = formula_libhomo, design = gss_all_filtered, 
##     family = quasibinomial())
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     -1.455e+02  3.521e+01  -4.133 4.40e-05 ***
## YEAR                             7.272e-02  1.747e-02   4.162 3.90e-05 ***
## RACEBlack                       -7.657e-01  1.113e-01  -6.882 2.45e-11 ***
## RACEOther                       -1.426e-01  1.493e-01  -0.955   0.3399    
## AGE                             -1.974e-02  2.873e-03  -6.868 2.66e-11 ***
## EDUC                             2.213e-01  1.513e-02  14.627  < 2e-16 ***
## SEXORNTBisexual                 -8.462e-01  6.324e-01  -1.338   0.1816    
## SEXORNTHeterosexual or Straight -1.074e+00  5.224e-01  -2.057   0.0404 *  
## Veteran_StatusVeteran            3.488e-02  1.578e-01   0.221   0.8252    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasibinomial family taken to be 1.007693)
## 
## Number of Fisher Scoring iterations: 5
library(survey)

# Expanded model including more independent variables
expanded_model_homosex <- svyglm(NUMERIC_HOMOSEX ~ YEAR + RACE + AGE + EDUC + SEXORNT + Veteran_Status, design = gss_all_filtered)

# Display the summary of the expanded model to see the coefficients, standard errors, and p-values
summary(expanded_model_homosex)
## 
## Call:
## svyglm(formula = NUMERIC_HOMOSEX ~ YEAR + RACE + AGE + EDUC + 
##     SEXORNT + Veteran_Status, design = gss_all_filtered)
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     -1.001e+02  1.847e+01  -5.418 1.07e-07 ***
## YEAR                             5.117e-02  9.175e-03   5.578 4.64e-08 ***
## RACEBlack                       -6.355e-01  6.114e-02 -10.393  < 2e-16 ***
## RACEOther                       -1.740e-01  7.316e-02  -2.379   0.0179 *  
## AGE                             -1.579e-02  1.447e-03 -10.915  < 2e-16 ***
## EDUC                             1.055e-01  6.965e-03  15.149  < 2e-16 ***
## SEXORNTBisexual                 -1.613e-01  1.091e-01  -1.478   0.1402    
## SEXORNTHeterosexual or Straight -9.050e-01  6.923e-02 -13.073  < 2e-16 ***
## Veteran_StatusVeteran           -1.749e-01  7.916e-02  -2.209   0.0278 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.732818)
## 
## Number of Fisher Scoring iterations: 2
# Assuming 'gss_all_filtered' is your survey design object created using svydesign

# Ensure MARHOMO is an ordered factor
gss_all_filtered$data$MARHOMO <- factor(gss_all_filtered$data$MARHOMO,
                                        levels = c("Strongly agree", "Agree", "Neither agree nor disagree", "Disagree", "Strongly disagree"),
                                        ordered = TRUE)

# Convert the ordered factor MARHOMO to numeric
gss_all_filtered$data$NUMERIC_MARHOMO <- as.numeric(gss_all_filtered$data$MARHOMO)

# Update the survey design object to include the new numeric variable
gss_all_filtered <- update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
# Expanded model with NUMERIC_MARHOMO as the dependent variable
expanded_model_marhomo <- svyglm(NUMERIC_MARHOMO ~ YEAR + RACE + AGE + EDUC + SEXORNT + Veteran_Status, design = gss_all_filtered)

# Display the summary of the expanded model
summary(expanded_model_marhomo)
## 
## Call:
## svyglm(formula = NUMERIC_MARHOMO ~ YEAR + RACE + AGE + EDUC + 
##     SEXORNT + Veteran_Status, design = gss_all_filtered)
## 
## Survey design:
## update(gss_all_filtered, NUMERIC_MARHOMO = as.numeric(MARHOMO))
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     188.221012  16.023382  11.747  < 2e-16 ***
## YEAR                             -0.092513   0.007955 -11.629  < 2e-16 ***
## RACEBlack                         0.395275   0.066630   5.932 6.71e-09 ***
## RACEOther                         0.067071   0.071931   0.932   0.3517    
## AGE                               0.019887   0.001384  14.374  < 2e-16 ***
## EDUC                             -0.093066   0.006695 -13.900  < 2e-16 ***
## SEXORNTBisexual                   0.154938   0.118804   1.304   0.1930    
## SEXORNTHeterosexual or Straight   0.996683   0.089983  11.076  < 2e-16 ***
## Veteran_StatusVeteran             0.174378   0.080855   2.157   0.0317 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 1.919018)
## 
## Number of Fisher Scoring iterations: 2