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