#Creating sampling weight
#Opening CSV data file
library(epicalc)
setwd((redacted))
use("Suppl1_data.csv")
des()
   
region_pop <- ifelse(region=="1 - Bangkok and Metro", 7375247,
ifelse(region=="2 - Central and Eastern", 12188254,
ifelse(region=="3 - Northern", 4793608,
ifelse(region=="4 - Western", 2472978,
ifelse(region=="5 - Northeastern", 16965035,
ifelse(region=="6 - Southern", 6930355, NA))))))
label.var(region_pop, "")

region_n <- ifelse(region=="1 - Bangkok and Metro", 323,
ifelse(region=="2 - Central and Eastern", 524,
ifelse(region=="3 - Northern", 208,
ifelse(region=="4 - Western", 110,
ifelse(region=="5 - Northeastern", 725,
ifelse(region=="6 - Southern", 301, NA))))))

swpop <- region_pop/region_n
label.var(swpop,"")
 
#Table 1. General Characteristics of the Study Participants
#Gender
#Male
#Female
sex <- ifelse(a1==1, "1 - Male", ifelse(a1==2, "2 - Female", ifelse(a1==3, "3 - Others", NA)))
label.var(sex, "")

sex_female <- ifelse(sex=="2 - Female", 1, 0)
label.var(sex_female, "")

sex_others <- ifelse(sex=="3 - Others", 1, 0)
label.var(sex_others, "")



#Age in years (mean ± SE)
age <- a3
age <- ifelse(age==88 | age == 99 | age > 120, NA, age)
label.var(age, "")

#Marital Status
#Single
#Married with child(ren)
#Married, no children
#Widowed/divorced/separated
marital_status <- ifelse(a4==1, "1 - Single", 
ifelse(a4==2, "2 - Married",
ifelse(a4==3, "3 - Cohabitation",
ifelse(a4==4, "4 - Divorced or Widowed",
ifelse(a4==5, "5 - It's complicated", NA)))))
label.var(marital_status, "")
tab1(marital_status)

marital2 <- ifelse(marital_status=="2 - Married", 1, 0)
marital3 <- ifelse(marital_status=="3 - Cohabitation", 1, 0)
marital4 <- ifelse(marital_status=="4 - Divorced or Widowed", 1, 0)
marital5 <- ifelse(marital_status=="5 - It's complicated", 1, 0)
label.var(marital2, "")
label.var(marital3, "")
label.var(marital4, "")
label.var(marital5, "")


#Highest Level of Education Completed
#Junior High School or Lower
#High School or Equivalent
#Some college (studied but no degree)
#Associate's degree
#Bachelor's Degree
#Graduate Degree
edu <- a5
edu_recat <- ifelse(edu <= 1, "1 - Junior HS or lower", 
ifelse(edu==2 | edu==3, "2 - HS or equivalent",
ifelse(edu==5, "3 - Some college or currently studying",
ifelse(edu==4, "4 - Assoc degree",
ifelse(edu==6, "5 - Bachelors degree",
ifelse(edu==7 | edu==8, "6 - Graduate degree", NA))))))
tab1(edu_recat)
label.var(edu_recat,"")

   #dummy variables (junior HS or lower as reference)
   #edu2 (high school or equivalent)
edu2 <- ifelse(edu_recat=="2 - HS or equivalent",1,0)
   #edu34 (some college but less than bachelor's degree)
edu34 <- ifelse(edu_recat=="3 - Some college or currently studying" | edu_recat=="4 - Assoc degree", 1, 0)
   #edu56 (bachelor's degree or higher)
edu56 <- ifelse(edu_recat=="5 - Bachelors degree" | edu_recat=="6 - Graduate degree", 1, 0)
label.var(edu2, "")
label.var(edu34, "")
label.var(edu56, "")

#Occupation
occ <- a6
occ <- ifelse(occ==10 | occ==11 | occ==16, "Cat 1",
ifelse(occ==3 | occ==5 | occ==7 | occ==9, "Cat 2", 
ifelse(occ==1 | occ==2 | occ==4 | occ==8 | occ==12, "Cat 3", NA)))
tab1(occ)
label.var(occ,"")
  #dummy variables
occ2 <- ifelse(occ=="Cat 2", 1, 0); label.var(occ2,"")
occ3 <- ifelse(occ=="Cat 3", 1, 0); label.var(occ3,"")




#Personal Monthly Income
#No more than 5,000 THB
#5,001 to 9,999 THB
#10,000 to 14,999 THB
#15,000 to 19,999 THB
#20,000 to 24,999 THB
#25,000 to 29,999 THB
#30,000 to 34,999 THB
#35,000 to 39,999 THB
#40,000 to 44,999 THB
#45,000 to 49,999 THB
#50,000 THB or more
income_recat <- a7
income_recat <- ifelse(income_recat==1, "01 - No more than 5k",
ifelse(income_recat==2, "02 - More than 5k to less than 10k",
ifelse(income_recat==3, "03 - 10k to less than 15k",
ifelse(income_recat==4, "04 - 15k to less than 20k",
ifelse(income_recat==5, "05 - 20k to less than 25k",
ifelse(income_recat==6, "06 - 25k to less than 30k",
ifelse(income_recat==7, "07 - 30k to less than 35k",
ifelse(income_recat==8, "08 - 35k to less than 40k",
ifelse(income_recat==9, "09 - 40k to less than 45k",
ifelse(income_recat==10, "10 - 45k to less than 50k",
ifelse(income_recat==11, "11 - 50k or more", NA)))))))))))
label.var(income_recat,"")
tab1(income_recat)


   #dummy (following the 30-30-30-10 ratio, arbitrarily)
   #Bottom 30% (income categories 1 and 2, making less than 10k per month)
income12 <- ifelse(income_recat=="01 - No more than 5k" | income_recat=="02 - More than 5k to less than 10k", 1, 0)

   #The next 30% (income category 3, making 10k to less than 15k per month)
income3 <- ifelse(income_recat=="03 - 10k to less than 15k", 1, 0)

   #The next 30% (income categories 4 and 5, making 15k to less than 25k per month)
income45 <- ifelse(income_recat=="04 - 15k to less than 20k" | income_recat=="05 - 20k to less than 25k", 1, 0)

   #Top 10% (income categories 6 thru 11, making 25k or more per month)
income611<- ifelse(income_recat=="06 - 25k to less than 30k" | income_recat=="07 - 30k to less than 35k" | 
income_recat=="08 - 35k to less than 40k" | income_recat=="09 - 40k to less than 45k" | 
income_recat=="10 - 45k to less than 50k" | income_recat=="11 - 50k or more", 1, 0)

label.var(income12, "")
label.var(income3, "")
label.var(income45, "")
label.var(income611, "")


#Religion
religion <- a9
religion <- ifelse(religion==1, "1 - Islam",
ifelse(religion==2, "2 - Buddhism", 
ifelse(religion==3, "3 - Christianity",
ifelse(religion==4, "4 - Others", NA))))
label.var(religion, "")
   #dummy variables (buddhists vs. non-buddhists)
non_buddhist <- ifelse(religion=="2 - Buddhism", 0, 1)
label.var(non_buddhist,"")



#Smoking status (not including electronic cigarettes)
#Never smokers (never did or smoked less than 100 cigs)
#Former smokers (smoked more than 100 cigs but not in past 12 months)
#Current smokers (smokedmore than 100 cigs in past 12 months)
ever_smoke <- b6_1
ever_smoke <- ifelse(ever_smoke==3, 1, 0)
  
sm_status <- ifelse(b6_2==888, "1 - Never smokers", "Unk")
sm_status <- ifelse(b6_1==2, "1 - Never smokers", sm_status)
sm_status <- ifelse(b6_1==3 & b6_2==1, "2 - Former smokers", sm_status)
sm_status <- ifelse(b6_1==3 & b6_2==2, "3 - Current smokers", sm_status)
sm_status <- ifelse(b6_1==3 & b6_2==3, "3 - Current smokers", sm_status)
sm_status <- ifelse(b6_1==3 & b6_2==4, "3 - Current smokers", sm_status)
sm_status <- ifelse(b6_2==999, NA, sm_status)
tab1(sm_status)
label.var(sm_status, "")
   #dummy variables
sm_former <- ifelse(sm_status=="2 - Former smokers", 1, 0)
sm_current <- ifelse(sm_status=="3 - Current smokers", 1, 0)
label.var(sm_former, "")
label.var(sm_current, "")



#Cannabis Use Status
#Never Users
#Former Users (used but not in past 12 months)
#Current Users (used within past 12 months)
   #First, coding for never vs. ever users
ever_users <- b1_1
   #Replacing refuse to answer as NA
ever_users <- ifelse(ever_users==99, NA, ever_users)
   #Coding proper for never vs. ever users
ever_users <- ifelse(ever_users==0, "0 - Never users", "1 - Ever users")
tab1(ever_users)

   #Second, coding for the most recent occasion
last_use <- b1_5
   #replacing missing values (888) with the never users and replacing the others with NA
last_use <- ifelse(b1_1==0, 0, last_use)
last_use <- ifelse(last_use==888, NA, last_use)
   #replacing "don't know" and "refuse to answer" as NA
last_use <- ifelse(last_use==88, NA, last_use)
last_use <- ifelse(last_use==99, NA, last_use)
tab1(last_use)

   #Coding for cannabis status based on the existing responses in the variable 'last_use'
cannabis_status <- ifelse(last_use==0, "0 - Never users",
ifelse(last_use==1, "1 - Former users", "2 - Current users"))
tab1(cannabis_status)
label.var(cannabis_status, "")

 
   #Coding as dummy variables for regression analyses
   #former vs. never users (current users treated as NA)
cannabis_status_former <- ifelse(cannabis_status=="1 - Former users", 1, 
ifelse(cannabis_status=="0 - Never users", 0, NA))
label.var(cannabis_status_former, "")

   #current vs. never users (former users treated as NA)
cannabis_status_current <- ifelse(cannabis_status=="2 - Current users", 1, 
ifelse(cannabis_status=="0 - Never users", 0, NA))
label.var(cannabis_status_current, "")

 
#Age of initiation
age_initiation <- b1_3_year
   #Excluding those who answered "Don't know" and refused to answer from the analyses
age_initiation <- ifelse(age_initiation>70, NA, age_initiation)
tab1(age_initiation)

#As a considerable portion (67 of 328, or 20.4%) of participants who were current cannabis users did not indicate their age of cannabis initiation, and cannabis initiation is one of the key determinants for cannabis use disorder and cannabis dosage...
#We decided to impute the median age of initiation (25 years of age) as the proxy age of initiation. 
#For participants who were younger than the imputed age of initiation, we imputed the age of 20 years as the proxy age of initiation. 
age_initiation[is.na(age_initiation)] <- 25
age_initiation <- ifelse(age<age_initiation, 20, age_initiation)
plot(age, age_initiation)
   #Participants who were never users never initiated, so we convert the age back to NA
age_initiation <- ifelse(cannabis_status=="0 - Never users", NA, age_initiation)
summary(age_initiation)
label.var(age_initiation, "")

   #Age of initiation without imputation
age_initiation_unimputed <- b1_3_year
   #Excluding those who answered "Don't know" and refused to answer from the analyses
age_initiation_unimputed <- ifelse(age_initiation_unimputed >70, NA, age_initiation_unimputed)
tab1(age_initiation_unimputed)
   #Participants who were never users never initiated, so we convert the age back to NA
age_initiation_unimputed <- ifelse(cannabis_status=="0 - Never users", NA, age_initiation_unimputed)
summary(age_initiation_unimputed)
label.var(age_initiation_unimputed , "")



#Past year initiation (Initiated within the past year (age of initiation differed from current age by no more than 1 year))
yearsago_initiate <- age-age_initiation
pastyear_initiate <- ifelse(yearsago_initiate<=1, 1, 0)
tab1(pastyear_initiate)
label.var(pastyear_initiate, "")



#Main purpose of cannabis use (among former and current users only)
#(if b1_4==99, recode as NA)
#Recreation (for enjoyment / relaxation / socialization) (b1_4==1)
#Other non-medical purposes (appetite, sleep, others such as cooking or experimentation) (b1_4==3 |b1_4==4 | b1_4==5)
#Medical purpose (pain relief, cancer, seizures) (b1_4==2)
cannabis_purpose <- b1_4
cannabis_purpose <- ifelse(cannabis_purpose>=99, NA, cannabis_purpose)
cannabis_purpose <- ifelse(cannabis_purpose==1, "1 - Recreation",
ifelse(cannabis_purpose==3 | cannabis_purpose==4 | cannabis_purpose==5, "2 - Other non-med", "3 - Medical"))
tab1(cannabis_purpose)
label.var(cannabis_purpose, "")


#Method of cannabis use (if b1_8==99, recode as NA)
cannabis_method <- b1_8
cannabis_method <- ifelse(cannabis_method==99, NA, cannabis_method)

#Ingestion (eaten or drank) (b1_8==5)
#Mixed with cigarette (b1_8==2)
#Rolled in a joint (b1_8==1)
#Smoked with water pipes (b1_8==4)
#Smoked in pipes or dried bong (b1_8==3)
#Others (infusion drop, capsules, spray, inhalants) (b1_8==6; b1_8==7; b1_8==8; b1_8==9; b1_8==10; b1_8==11)
cannabis_method <- ifelse(cannabis_method==5, "1 - Ingestion",
ifelse(cannabis_method==2, "2 - Mixed with cigs",
ifelse(cannabis_method==1, "3 - Rolled in joint",
ifelse(cannabis_method==4, "4 - Smoked water pipes",
ifelse(cannabis_method==3, "5 - Smoked in dry bongs",
ifelse(cannabis_method==6 | cannabis_method==7 | cannabis_method==8 | cannabis_method==9 | cannabis_method==10 | cannabis_method==11, "6 - Others", NA))))))
tab1(cannabis_method)
label.var(cannabis_method, "")
 
#Density of cannabis outlets within a given distance from place of residence
#Density within 400 meters 
#None
#1-2 outlets
#3 or more outlets
density_400m <- c1_6_400m
density_400m <- ifelse(density_400m==99, NA, density_400m)
density_400m <- ifelse(density_400m==88, 0, density_400m)
tab1(density_400m)

   #Categorized
density_400m_cats <- ifelse(density_400m==0, "0 - None",
ifelse(density_400m==1, "1 - One",
ifelse(density_400m==2, "2 - Two",
ifelse(density_400m>=3, "3 - Three or more", NA))))
tab1(density_400m_cats)
label.var(density_400m_cats, "")

   #Dummy variables (with "none" as reference group)
density_400m_1 <- ifelse(density_400m_cats=="1 - One", 1, 0)
density_400m_2 <- ifelse(density_400m_cats=="2 - Two", 1, 0)
density_400m_3 <- ifelse(density_400m_cats=="3 - Three or more", 1, 0)
label.var(density_400m_1, "")
label.var(density_400m_2, "")
label.var(density_400m_3, "")




#Density within 800 meters 
#None
#1-2 outlets
#3 or more outlets
density_800m <- c1_7_800m
density_800m <- ifelse(density_800m==99, NA, density_800m)
density_800m <- ifelse(density_800m==88, 0, density_800m)
tab1(density_800m)

   #Categorized
density_800m_cats <- ifelse(density_800m==0, "0 - None",
ifelse(density_800m==1, "1 - One",
ifelse(density_800m==2, "2 - Two",
ifelse(density_800m>=3, "3 - Three or more", NA))))
tab1(density_800m_cats)
label.var(density_800m_cats, "")

   #Dummy variables (with "none" as reference group)
density_800m_1 <- ifelse(density_800m_cats=="1 - One", 1, 0)
density_800m_2 <- ifelse(density_800m_cats=="2 - Two", 1, 0)
density_800m_3 <- ifelse(density_800m_cats=="3 - Three or more", 1, 0)
label.var(density_800m_1, "")
label.var(density_800m_2, "")
label.var(density_800m_3, "")


#Density within 1200 meters 
#None
#1-2 outlets
#3 or more outlets
density_1200m <- c1_8_1200m
density_1200m <- ifelse(density_1200m==99, NA, density_1200m)
density_1200m <- ifelse(density_1200m==88, 0, density_1200m)
tab1(density_1200m)

   #Categorized
density_1200m_cats <- ifelse(density_1200m==0, "0 - None",
ifelse(density_1200m==1, "1 - One",
ifelse(density_1200m==2, "2 - Two",
ifelse(density_1200m>=3, "3 - Three or more", NA))))
tab1(density_1200m_cats)
label.var(density_1200m_cats, "")

   #Dummy variables (with "none" as reference group)
density_1200m_1 <- ifelse(density_1200m_cats=="1 - One", 1, 0)
density_1200m_2 <- ifelse(density_1200m_cats=="2 - Two", 1, 0)
density_1200m_3 <- ifelse(density_1200m_cats=="3 - Three or more", 1, 0)
label.var(density_1200m_1, "")
label.var(density_1200m_2, "")
label.var(density_1200m_3, "")

 
#Density within 1600 meters 
#None
#1-2 outlets
#3 or more outlets
density_1600m <- c1_9_1600m
density_1600m <- ifelse(density_1600m==99, NA, density_1600m)
density_1600m <- ifelse(density_1600m==88, 0, density_1600m)
tab1(density_1600m)

   #Categorized
density_1600m_cats <- ifelse(density_1600m==0, "0 - None",
ifelse(density_1600m==1, "1 - One",
ifelse(density_1600m==2, "2 - Two",
ifelse(density_1600m>=3, "3 - Three or more", NA))))
tab1(density_1600m_cats)
label.var(density_1600m_cats, "")

   #Dummy variables (with "none" as reference group)
density_1600m_1 <- ifelse(density_1600m_cats=="1 - One", 1, 0)
density_1600m_2 <- ifelse(density_1600m_cats=="2 - Two", 1, 0)
density_1600m_3 <- ifelse(density_1600m_cats=="3 - Three or more", 1, 0)
label.var(density_1600m_1, "")
label.var(density_1600m_2, "")
label.var(density_1600m_3, "")





#Cannabis Use Disorder Status (among current users only)
   #Treat non-responses as NA. 
   #First, create the score variables
cud1 <- b3_1
cud2 <- b3_2
cud3 <- b3_3
cud4 <- b3_4
cud5 <- b3_5
cud6 <- b3_6
cud7 <- b3_7
cud8 <- b3_8
 
  #Second, Treat skips of non-current users (those who answered no use of cannabis in past 12 months in Question 1) as NA
cud1 <- ifelse(cud1==0, NA, cud1)
cud2 <- ifelse(cud2==888, NA, cud2)
cud3 <- ifelse(cud3==888, NA, cud3)
cud4 <- ifelse(cud4==888, NA, cud4)
cud5 <- ifelse(cud5==888, NA, cud5)
cud6 <- ifelse(cud6==888, NA, cud6)
cud7 <- ifelse(cud7==888, NA, cud7)
cud8 <- ifelse(cud8==888, NA, cud8)

   #Third, treat non-responses as NA
cud1 <- ifelse(cud1==9, NA, cud1)
cud2 <- ifelse(cud2==9, NA, cud2)
     #for CUD Question 2, also recoded "No symptoms of being high" as 0 (less than 1 hour)
cud2 <- ifelse(cud2==10, 0, cud2)
cud3 <- ifelse(cud3==9, NA, cud3)
cud4 <- ifelse(cud4==9, NA, cud4)
cud5 <- ifelse(cud5==9, NA, cud5)
cud6 <- ifelse(cud6==9, NA, cud6)
cud7 <- ifelse(cud7==9, NA, cud7)
cud8 <- ifelse(cud8==9, NA, cud8)

   #Fourth, calculate the score proper based on the original CUDIT-R instrument
     #Item 1: Frequency of cannabis use (CUD1): 
          #Never =0 point; 
          #Monthly or less = 1 point;
          #2-4 times a month = 2 points; 
          #2-3 times a week = 3 points; 
          #4 or more times per week = 4 points
cud1_score <- ifelse(cud1==0, 0, ifelse(cud1>=1 & cud1<=4, 1, ifelse(cud1==5, 2, ifelse(cud1==6, 3, 4))))
tab1(cud1_score)
label.var(cud1_score, "")

     #Item 2 thru 7 are already appropriate
cud2_score <- cud2
cud3_score <- cud3
cud4_score <- cud4
cud5_score <- cud5
cud6_score <- cud6
cud7_score <- cud7
label.var(cud2_score, "")
label.var(cud3_score, "")
label.var(cud4_score, "")
label.var(cud5_score, "")
label.var(cud6_score, "")
label.var(cud7_score, "")

     #Item 8: Thought of cutting down or stopping cannabis use (this item is also already appropriate
cud8_score <- cud8
label.var(cud8_score, "")

  #Third, calculate the score (variable name 'cud_core')
cud_score <- cud1_score + cud2_score + cud3_score + cud4_score + cud5_score + cud6_score + cud7_score + cud8_score
tab1(cud_score)
label.var(cud_score,"")


#Fourth, categorize according to the score
#คณะผู้วิจัยจะใช้แบบสอบถาม Cannabis Use Disorder Identification Test - Revised (CUDIT-R) 
#ซึ่งเป็นแบบสอบถามที่ใช้คัดกรองอาการติดกัญชา มีแบบสอบถาม 8 ข้อ คะแนนตั้งแต่ 0-24 คะแนน โดยมีจุดตัดของ hazardous use ที่ 8-11 คะแนน และมีจุดตัดของ possible cannabis use disorder ที่ 12-24 คะแนน 
#Non-hazardous use (1-7 points)
#Hazardous use (8-11 points)
#Possible cannabis use disorder (12 points or higher)
#Cannabis use amount within past year (grams per person per day of use, current users only)
cud_status <- ifelse(cud_score <= 7, "1 - Non-hazardous use",
ifelse(cud_score >=8 & cud_score <= 11, "2 - Hazardous use", "3 - Possible CUD"))
tab1(cud_status)
label.var(cud_status, "")

   #Create binary variable for cross-tabulation and regression analyses
cud_haz_disorder <- ifelse(cud_status=="2 - Hazardous use" | cud_status=="3 - Possible CUD", 1, 0)
label.var(cud_haz_disorder, "")


#Amount of cannabis use per day of use
#Amount of cannabis user per day of use was calculated based on the response to the questions: 1) "B1_12a.  In each cannabis use session, in general, how many grams cannabis do you use yourself?" and; 2) "B1_12b. On the day that you use cannabis, in general, how many times do you use cannabis per day?". After we treated the "don't know" and "refuse to answer" of each question as NA, we multipled the two numeric responses in order to obtain the approximately amount of cannabis use by the participant (in grams per day of use). 

cannabis_use_grams <- b1_12a_gram
cannabis_useday_times <- b1_12b_times
cannabis_use_grams <- ifelse(cannabis_use_grams==888 | cannabis_use_grams==999, NA, cannabis_use_grams)
cannabis_useday_times <- ifelse(cannabis_useday_times==888 | cannabis_useday_times==999, NA, cannabis_useday_times)
label.var(cannabis_use_grams,"")
label.var(cannabis_useday_times,"")

cannabis_useday_grams <- cannabis_use_grams * cannabis_useday_times
label.var(cannabis_useday_grams, "")

#Limiting the analyses only to current users
cannabis_useday_grams_current <- ifelse(cannabis_status=="2 - Current users", cannabis_useday_grams, NA)
label.var(cannabis_useday_grams_current, "")

#Table 4. Cannabis Use Amount (Grams per Person per Day of Use) within past year among Participants by Density of Cannabis Outlets (Among Current Users Only)
summ(cannabis_useday_grams_current, by=density_400m_cats)
summ(cannabis_useday_grams_current, by=density_800m_cats)
summ(cannabis_useday_grams_current, by=density_1200m_cats)
summ(cannabis_useday_grams_current, by=density_1600m_cats)


#SUPPL Table 1
#Supplementary Table 1. Cannabis use status of study participants who reported vs. did not report the number of outlets (answered “don’t know”) within each given distance (row weighted percent ± standard errors, unless otherwise indicated)

#Creating variable density_dontknow_400m for those who reported the number of outlets within 400m radius (0) vs. answered "don't know" (1), excluding those who refused to answer from the analyses
   #Repeating the same codes for the other distances
density_dontknow_400m <-   c1_6_400m
density_dontknow_800m <-   c1_7_800m
density_dontknow_1200m <- c1_8_1200m
density_dontknow_1600m <- c1_9_1600m

density_dontknow_400m <- ifelse(density_dontknow_400m==99, NA, density_dontknow_400m)
density_dontknow_800m <- ifelse(density_dontknow_800m==99, NA, density_dontknow_800m)
density_dontknow_1200m <- ifelse(density_dontknow_1200m==99, NA, density_dontknow_1200m)
density_dontknow_1600m <- ifelse(density_dontknow_1600m==99, NA, density_dontknow_1600m)

density_dontknow_400m <- ifelse(density_dontknow_400m==88, 1, 0)
density_dontknow_800m <- ifelse(density_dontknow_800m==88, 1, 0) 
density_dontknow_1200m <- ifelse(density_dontknow_1200m==88, 1, 0) 
density_dontknow_1600m <- ifelse(density_dontknow_1600m==88, 1, 0) 

tab1(density_dontknow_400m)
tab1(density_dontknow_800m)
tab1(density_dontknow_1200m)
tab1(density_dontknow_1600m)

label.var(density_dontknow_400m, "")
label.var(density_dontknow_800m, "")
label.var(density_dontknow_1200m, "")
label.var(density_dontknow_1600m, "")
 
tabpct(density_dontknow_400m, cannabis_status)
tabpct(density_dontknow_800m, cannabis_status)
tabpct(density_dontknow_1200m, cannabis_status)
tabpct(density_dontknow_1600m, cannabis_status)


#SUPPL Table 2 (Sensitivity Analysis)
#Supplementary Table 2. Density (without imputation of "don't know" as 0 vs. Cannabis use status)
sen2_400m 	<- ifelse(density_dontknow_400m==0, density_400m_cats, NA)
sen2_400m_1 	<- ifelse(density_dontknow_400m==0, density_400m_1, NA)
sen2_400m_2 	<- ifelse(density_dontknow_400m==0, density_400m_2, NA)
sen2_400m_3 	<- ifelse(density_dontknow_400m==0, density_400m_3, NA)

sen2_800m 	<- ifelse(density_dontknow_800m==0, density_800m_cats, NA)
sen2_800m_1 	<- ifelse(density_dontknow_800m==0, density_800m_1, NA)
sen2_800m_2 	<- ifelse(density_dontknow_800m==0, density_800m_2, NA)
sen2_800m_3 	<- ifelse(density_dontknow_800m==0, density_800m_3, NA)

sen2_1200m 	<- ifelse(density_dontknow_1200m==0, density_1200m_cats, NA)
sen2_1200m_1 	<- ifelse(density_dontknow_1200m==0, density_1200m_1, NA)
sen2_1200m_2 	<- ifelse(density_dontknow_1200m==0, density_1200m_2, NA)
sen2_1200m_3 	<- ifelse(density_dontknow_1200m==0, density_1200m_3, NA)

sen2_1600m 	<- ifelse(density_dontknow_1600m==0, density_1600m_cats, NA)
sen2_1600m_1 	<- ifelse(density_dontknow_1600m==0, density_1600m_1, NA)
sen2_1600m_2 	<- ifelse(density_dontknow_1600m==0, density_1600m_2, NA)
sen2_1600m_3 	<- ifelse(density_dontknow_1600m==0, density_1600m_3, NA)

tabpct(sen2_400m, cannabis_status)
tabpct(sen2_800m, cannabis_status)
tabpct(sen2_1200m, cannabis_status)
tabpct(sen2_1600m, cannabis_status)
 

label.var(sen2_400m, "")
label.var(sen2_400m_1, "")
label.var(sen2_400m_2, "")
label.var(sen2_400m_3, "")
label.var(sen2_800m, "")
label.var(sen2_800m_1, "")
label.var(sen2_800m_2, "")
label.var(sen2_800m_3, "")
label.var(sen2_1200m, "")
label.var(sen2_1200m_1, "")
label.var(sen2_1200m_2, "")
label.var(sen2_1200m_3, "")
label.var(sen2_1600m, "")
label.var(sen2_1600m_1, "")
label.var(sen2_1600m_2, "")
label.var(sen2_1600m_3, "")

 
################################################
#Weighted analyses 20230606
library(survey)
design <- svydesign(id=~id, strata=~region, weights=~swpop, data=.data)



#TABLE 1
#Sex
svymean(~sex, design=design, data=.data, na.rm=TRUE)
  
#Age in years (mean ± SE)
svymean(~age, design=design, data=.data, na.rm=TRUE)

#Marital Status
#Single
#Married with child(ren)
#Married, no children
#Widowed/divorced/separated
svymean(~marital_status, design=design, data=.data, na.rm=TRUE)

#Highest Level of Education Completed
#Junior High School or Lower
#High School or Equivalent
#Some college (did not graduate) or currently studying
#Associate’s degree
#Bachelor's Degree
#Graduate Degree
svymean(~edu_recat, design=design, data=.data, na.rm=TRUE)

#Occupation categories
svymean(~occ, design=design, data=.data, na.rm=TRUE)

Personal Monthly Income
#No more than 5,000 THB
#5,001 to 9,999 THB
#10,000 to 14,999 THB
#15,000 to 19,999 THB
#20,000 to 24,999 THB
#25,000 to 29,999 THB
#30,000 to 34,999 THB
#35,000 to 39,999 THB
#40,000 to 44,999 THB
#45,000 to 49,999 THB
#50,000 THB or more
svymean(~income_recat, design=design, data=.data, na.rm=TRUE)

#Religion
svymean(~religion, design=design, data=.data, na.rm=TRUE)

#Smoking status
svymean(~sm_status, design=design, data=.data, na.rm=TRUE)
 
 
#Density of cannabis outlets within a given distance from place of residence*
#400 meters (walking for 5 minutes)
svymean(~density_400m_cats, design=design, data=.data, na.rm=TRUE)

#800 meters (walking for 10 minutes)
svymean(~density_800m_cats, design=design, data=.data, na.rm=TRUE)
 
#1200 meters (walking for 15 minutes)
svymean(~density_1200m_cats, design=design, data=.data, na.rm=TRUE)
 
#1600 meters
svymean(~density_1600m_cats, design=design, data=.data, na.rm=TRUE)
 
#Cannabis Use Status
#Never Users
#Former Users (used but not in past 12 months)
#Current Users (used within past 12 months)
svymean(~cannabis_status, design=design, data=.data, na.rm=TRUE)
 
#Age of initiation (years) for former and current users (mean ± SE)
svymean(~age_initiation, design=design, data=.data, na.rm=TRUE)

#Initiated within past year
svymean(~pastyear_initiate, design=design, data=.data, na.rm=TRUE)

#Main purpose of cannabis use
svymean(~cannabis_purpose, design=design, data=.data, na.rm=TRUE)
 
#Method of cannabis use
svymean(~cannabis_method, design=design, data=.data, na.rm=TRUE)

 #Cannabis Use Disorder Status (among current users only)
#Non-hazardous use (0-7 points)
#Hazardous use (8-11 points)
#Possible cannabis use disorder (12 points or higher)
svymean(~cud_status, design=design, data=.data, na.rm=TRUE)



#Cannabis use amount within past year (grams per person per day of use, current users only)
svymean(~cannabis_useday_grams_current, design=design, data=.data, na.rm=TRUE)
 
 
#TABLE 2
   #Cross-tabulation with unweighted data first, to get the sample size for each row
#tabpct(density_400m_cats, cannabis_status)
#tabpct(density_800m_cats, cannabis_status)
#tabpct(density_1200m_cats, cannabis_status)
#tabpct(density_1600m_cats, cannabis_status)

#Distances vs. cannabis use status (row percents only)
svyby(~cannabis_status, ~density_400m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~density_800m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~density_1200m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~density_1600m_cats, design=design, svymean, na.rm=TRUE) 

#Crude Odds ratios with 95% CI
#At 400m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ density_400m_1 + density_400m_2 + density_400m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ density_400m_1 + density_400m_2 + density_400m_3, family=binomial, design = design))
summary(model2)

   #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ density_400m_1 + density_400m_2 + density_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
  
   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ density_400m_1 + density_400m_2 + density_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
 
#At 800m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ density_800m_1 + density_800m_2 + density_800m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ density_800m_1 + density_800m_2 + density_800m_3, family=binomial, design = design))
summary(model2)

   #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ density_800m_1 + density_800m_2 + density_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
  
   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ density_800m_1 + density_800m_2 + density_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)

#At 1200 m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ density_1200m_1 + density_1200m_2 + density_1200m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ density_1200m_1 + density_1200m_2 + density_1200m_3, family=binomial, design = design))
summary(model2)
 
   #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ density_1200m_1 + density_1200m_2 + density_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)

   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ density_1200m_1 + density_1200m_2 + density_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)


#At 1600m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ density_1600m_1 + density_1600m_2 + density_1600m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ density_1600m_1 + density_1600m_2 + density_1600m_3, family=binomial, design = design))
summary(model2)

 #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ density_1600m_1 + density_1600m_2 + density_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
  
   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ density_1600m_1 + density_1600m_2 + density_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)


 
#TABLE 3
   #Cross-tabulation with unweighted data first, to get the sample size for each row
#Distances vs. cannabis use status (row percents only)
svyby(~cud_haz_disorder, ~density_400m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cud_haz_disorder, ~density_800m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cud_haz_disorder, ~density_1200m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cud_haz_disorder, ~density_1600m_cats, design=design, svymean, na.rm=TRUE) 

#Crude Odds ratios with 95% CI
#At 400m
   #Crude OR (95% CI) for Hazardous or Potential disorder vs. Non-Hazardous use
model1<-(svyglm(cud_haz_disorder ~ density_400m_1 + density_400m_2 + density_400m_3, family=binomial, design = design))
summary(model1)

#*Adjusted for the participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use
#Adjusted OR
model2<-(svyglm(cud_haz_disorder ~ density_400m_1 + density_400m_2 + density_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, family=binomial, design = design))
summary(model2)



 
#At 800m
   #Crude OR (95% CI) for Hazardous or Potential disorder vs. Non-Hazardous use
model1<-(svyglm(cud_haz_disorder ~ density_800m_1 + density_800m_2 + density_800m_3, family=binomial, design = design))
summary(model1)

#Adjusted OR
model2<-(svyglm(cud_haz_disorder ~ density_800m_1 + density_800m_2 + density_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, family=binomial, design = design))
summary(model2)


#At 1200 m
   #Crude OR (95% CI) for Hazardous or Potential disorder vs. Non-Hazardous use
model1<-(svyglm(cud_haz_disorder ~ density_1200m_1 + density_1200m_2 + density_1200m_3, family=binomial, design = design))
summary(model1)

#Adjusted OR
model2<-(svyglm(cud_haz_disorder ~ density_1200m_1 + density_1200m_2 + density_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, family=binomial, design = design))
summary(model2)


#At 1600m
   #Crude OR (95% CI) for Hazardous or Potential disorder vs. Non-Hazardous use
model1<-(svyglm(cud_haz_disorder ~ density_1600m_1 + density_1600m_2 + density_1600m_3, family=binomial, design = design))
summary(model1)
 
 #Adjusted OR
model2<-(svyglm(cud_haz_disorder ~ density_1600m_1 + density_1600m_2 + density_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, family=binomial, design = design))
summary(model2)
 
 
#TABLE 4 Cannabis Use Amount (Grams per Person per Day of Use) within past year among Participants by Density of Cannabis Outlets (Among Current Users Only) (n=84)
   #Cross-tabulation with unweighted data first, to get the sample size for each row
#summ(cannabis_useday_grams_current, by=density_400m_cats)
#summ(cannabis_useday_grams_current, by=density_800m_cats)
#summ(cannabis_useday_grams_current, by=density_1200m_cats)
#summ(cannabis_useday_grams_current, by=density_1600m_cats)

svyby(~cannabis_useday_grams_current, ~density_400m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_useday_grams_current, ~density_800m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_useday_grams_current, ~density_1200m_cats, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_useday_grams_current, ~density_1600m_cats, design=design, svymean, na.rm=TRUE) 


#Beta (Linear Regression Coefficient) with 95% CI
#At 400m
   #Beta (Linear Regression Coefficient) with 95% CI
model1<-(svyglm(cannabis_useday_grams_current ~ density_400m_1 + density_400m_2 + density_400m_3, design = design))
summary(model1)

model2<-(svyglm(cannabis_useday_grams_current ~ density_400m_1 + density_400m_2 + density_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model2)



 
#At 800m
   #Beta (Linear Regression Coefficient) with 95% CI
model1<-(svyglm(cannabis_useday_grams_current ~ density_800m_1 + density_800m_2 + density_800m_3,  design = design))
summary(model1)

model2<-(svyglm(cannabis_useday_grams_current ~ density_800m_1 + density_800m_2 + density_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model2)



 
#At 1200 m
   #Beta (Linear Regression Coefficient) with 95% CI
model1<-(svyglm(cannabis_useday_grams_current ~ density_1200m_1 + density_1200m_2 + density_1200m_3,  design = design))
summary(model1)

model2<-(svyglm(cannabis_useday_grams_current ~ density_1200m_1 + density_1200m_2 + density_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model2)



 
#At 1600m
   #Beta (Linear Regression Coefficient) with 95% CI
model1<-(svyglm(cannabis_useday_grams_current ~ density_1600m_1 + density_1600m_2 + density_1600m_3,  design = design))
summary(model1)
 

model2<-(svyglm(cannabis_useday_grams_current ~ density_1600m_1 + density_1600m_2 + density_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model2)
 







######################
#SUPPL Table 1
#Supplementary Table 1. Cannabis use status of study participants who reported vs. did not report the number of outlets (answered “don’t know”) within each given distance (row weighted percent ± standard errors, unless otherwise indicated)
svyby(~cannabis_status, ~density_dontknow_400m, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~density_dontknow_800m, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~density_dontknow_1200m, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~density_dontknow_1600m, design=design, svymean, na.rm=TRUE) 
svychisq(~cannabis_status+density_dontknow_400m, design)
svychisq(~cannabis_status+density_dontknow_800m, design)
svychisq(~cannabis_status+density_dontknow_1200m, design)
svychisq(~cannabis_status+density_dontknow_1600m, design)

###########################
#SUPPL Table 2
#Distances vs. cannabis use status (row percents only)
svyby(~cannabis_status, ~sen2_400m, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~sen2_800m, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~sen2_1200m, design=design, svymean, na.rm=TRUE) 
svyby(~cannabis_status, ~sen2_1600m, design=design, svymean, na.rm=TRUE) 

#Crude Odds ratios with 95% CI
#At 400m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ sen2_400m_1 + sen2_400m_2 + sen2_400m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ sen2_400m_1 + sen2_400m_2 + sen2_400m_3, family=binomial, design = design))
summary(model2)

   #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ sen2_400m_1 + sen2_400m_2 + sen2_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
  
   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ sen2_400m_1 + sen2_400m_2 + sen2_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
 
#At 800m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ sen2_800m_1 + sen2_800m_2 + sen2_800m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ sen2_800m_1 + sen2_800m_2 + sen2_800m_3, family=binomial, design = design))
summary(model2)

   #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ sen2_800m_1 + sen2_800m_2 + sen2_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
  
   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ sen2_800m_1 + sen2_800m_2 + sen2_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)

#At 1200 m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ sen2_1200m_1 + sen2_1200m_2 + sen2_1200m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ sen2_1200m_1 + sen2_1200m_2 + sen2_1200m_3, family=binomial, design = design))
summary(model2)
 
   #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ sen2_1200m_1 + sen2_1200m_2 + sen2_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)

   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ sen2_1200m_1 + sen2_1200m_2 + sen2_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)


#At 1600m
   #Crude OR (95% CI) for Former vs. Never Users
model1<-(svyglm(cannabis_status_former ~ sen2_1600m_1 + sen2_1600m_2 + sen2_1600m_3, family=binomial, design = design))
summary(model1)

   #Crude OR (95% CI) for Current vs. Never Users
model2<-(svyglm(cannabis_status_current ~ sen2_1600m_1 + sen2_1600m_2 + sen2_1600m_3, family=binomial, design = design))
summary(model2)

 #Adjusted OR (95% CI) for Former vs. Never users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_former ~ sen2_1600m_1 + sen2_1600m_2 + sen2_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
  
   #Adjusted OR (95% CI) for Current vs. Never Users (participant's sex, age, tobacco smoking status, marital status, income, religion, occupation, educational level, and age of onset of cannabis use)
model1<-(svyglm(cannabis_status_current ~ sen2_1600m_1 + sen2_1600m_2 + sen2_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + religion + occ2 + occ3 + edu2 + edu34 + edu56, family=binomial, design = design))
summary(model1)
 
 
###############
#Supplementary Table 3. Adjusted Beta (95% CI) for cannabis use amount (grams per person per day of use) within the past year by density of cannabis outlets (among current users who smoked the cannabis) without imputation of age of initiation

#MODEL 1 (with imputation for age at initiation)
#Adj Beta (Linear Regression Coefficient) with 95% CI
#At 400m
model1<-(svyglm(cannabis_useday_grams_current ~ sen2_400m_1 + sen2_400m_2 + sen2_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model1)
  
#At 800m
model1<-(svyglm(cannabis_useday_grams_current ~ sen2_800m_1 + sen2_800m_2 + sen2_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model1)

#At 1200 m
model1<-(svyglm(cannabis_useday_grams_current ~ sen2_1200m_1 + sen2_1200m_2 + sen2_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model1)

#At 1600m
model1<-(svyglm(cannabis_useday_grams_current ~ sen2_1600m_1 + sen2_1600m_2 + sen2_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation, design = design))
summary(model1)
  


#MODEL 2 (without imputation for age at initiation)
#Adj Beta (Linear Regression Coefficient) with 95% CI
#At 400m
model2<-(svyglm(cannabis_useday_grams_current ~ sen2_400m_1 + sen2_400m_2 + sen2_400m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation_unimputed, design = design))
summary(model2)
  
#At 800m
model2<-(svyglm(cannabis_useday_grams_current ~ sen2_800m_1 + sen2_800m_2 + sen2_800m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation_unimputed, design = design))
summary(model2)

#At 1200 m
model2<-(svyglm(cannabis_useday_grams_current ~ sen2_1200m_1 + sen2_1200m_2 + sen2_1200m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation_unimputed, design = design))
summary(model2)

#At 1600m
model2<-(svyglm(cannabis_useday_grams_current ~ sen2_1600m_1 + sen2_1600m_2 + sen2_1600m_3 + sex_female + sex_others + age + sm_former + sm_current + marital2 + marital3 + marital4 + marital5 + income12 + income3 + income45 + income611 + non_buddhist + occ2 + occ3 + edu2 + edu34 + edu56 + age_initiation_unimputed, design = design))
summary(model2)