Combining yearly data sets to perform svydesign analysis (CDC-BRFSS)

[mohddar]

I am trying combine data sets for the year 2016 till 2020, i loaded datasets (brrfss2016vars.rda ),
options(survey.lonely.psu = "adjust")

load("brfss2016vars.Rda")
load("brfss2017vars.Rda")
load("brfss2018vars.Rda")
load("brfss2020vars.Rda")

str(brfss2016vars)

Gender SEX1

Age agecat

Depression dep_disorder

E-cig: Never, Current, Former

Smk: Never, Current, Former

library(survey)
options(survey.lonely.psu="adjust")

prev <-svymean(~female,design=svy,na.rm=TRUE) #svymean works for variables coded as 0 or 1

confint(svymean(~female,design = svy,na.rm = TRUE))

svymean(~smkstat,design=svy,na.rm=TRUE)

svymean(~currentsmoker,design=svy,na.rm=TRUE)

svy16 <-svydesign(id=~1,strata=~X_STSTR,weights=~X_LLCPWT,data=brfss2016vars) # specify survey design for BRFSS 2016
svy17 <-svydesign(id=~1,strata=~X_STSTR,weights=~X_LLCPWT,data=brfss2017vars) # specify survey design for BRFSS 2017
svy18 <-svydesign(id=~1,strata=~X_STSTR,weights=~X_LLCPWT,data=brfss2018vars) # specify survey design for BRFSS 2018
svy20 <-svydesign(id=~1,strata=~X_STSTR,weights=~X_LLCPWT,data=brfss2020vars) # specify survey design for BRFSS 2020

Generate age group-specific prevalences and confidence intervals across adult population

getprevsbyage <- function(groupvar, subpop){

prevs = NULL
gender = c("males","females")
ages =c("18.24","25.34","35.44","45.54","55.64","65")
depstat = c("notdep","dep")
years = c(2016,2017,2018,2020) #brfss years with ecig data
svydesigns = list(svy16, svy17, svy18, svy20)

for (g in c(1:2)){ # male, female
for (a in c(1:4)){
for (d in c(0:1)){
for (s in years){

    subsvy <- subset(svydesigns[s],(female==g & agecat==a & dep==d))
    
    p <- svymean(~neversmoker, design=subsvy, na.rm=TRUE)
    q <- svymean(~currentsmoker, design=subsvy, na.rm=TRUE)
    r <- svymean(~formersmoker, design=subsvy, na.rm=TRUE)
    
    prevs <- rbind(prevs, cbind(round(p,3),round(confint(p)[1],3), round(confint(p)[2],3), 
                               gender[g],ages[a],depstat[d+1], "neversmoker",years[s]), 
    cbind(round(q,3),round(confint(q)[1],3), round(confint(q)[2],3), 
                                 gender[g],ages[a],depstat[d+1],"currentsmoker",years[s]),
    cbind(round(r,3),round(confint(r)[1],3), round(confint(r)[2],3), 
                                 gender[g],ages[a],depstat[d+1],"formersmoker",years[s]))
  }
}

}
}
colnames(prevs) <- c("prev","lowCI","highCI","gender","agecat","depstat","smkstat","year")
prevs

However i am getting error in "for" loop won't give me prevalence? could you please help to find what's wrong in the code.
Thank you