/* Biostat 200 Lab 1 This lab will introduce some basic commands to * inspect a dataset * do some basic data management * demonstrate use of do-file storing your commands, and * make a log file to capture the output. The do-files we provide for these labs use a bit more than the most basic Stata programming, for the benefit of students who want to learn it. For example, * loops, for performing essentially the same command repeatedly with minor changes to the inputs * local and global variables to store a number or a string, which can then be plugged in elsewhere as needed. If this is not of interest, you can simply run the code we provide, focusing on understanding the result. To use this do-file, we suggest that you highlight the commands, or blocks of commands, between each block of comments and instructions, then press the do button in the upper right corner of the do-file editor. In some cases, it is crucial to run a whole block of code containing a loop together for it to work, starting from the line with a foreach or forvalues command, and running through the closing bracket }, which is on a separate line. First, modify the next line so it gives the path to the directory where you have put this do-file and the dataset lab1.dta. */ cd "~/work/mb/atcr/200/labs" /* Keeping a log file is a good way to record the results of the do-file. We first make sure there is no existing log file open, then start a new one. The capture prefix prevents Stata from stopping if there is no log file open. This is needed when you are running a do-file repeatedly to debug it. */ capture log close log using "lab 1", replace /* Read in the data posted on the class web site. These data come from a study of voluntary counseling and testing (VCT) for HIV at Mulago Hospital in Kampala, Uganda. */ use vct.dta, clear *Stata lets us look at a spreadsheet of the data: */ browse /* First see what variables are in the dataset. Check how many observations and variables the dataset includes. */ describe /* Now get brief summary statistics for all the variables in the dataset. The full command we'll use is summarize, but we can shorten that to save typing. */ sum /* Now do a little data management to make the data easier to use. First assign some variable labels */ label var age "Age" label var auditc "Audit score, last 3 months" label var cd4 "CD4 count" label var cd4_cat "CD4 count" label var hiv "HIV status" label var hungry "How often household members go hungry" label var lastalc "Last time used alcohol" label var female "female sex" label var n_adults "# adults in household" label var n_child "# people <18yo in household" label var n_support "# people supported financially" label var religion "Religion" /* Now add some missing value labels. A easy-to-remember convention is to use the variable name as the name of the value label. (Stata doesn't have any problem with this). Other easy-to-remember conventions will also work (e.g., add an underscore to the variable name).*/ label define yesno 0 "no" 1 "yes" label values female yesno label define sex 1 "men" 2 "women" label values sex sex label define hiv 0 "HIV-" 1 "HIV+" label values hiv hiv label define cd4_cat 0 "<50" 1 "51-250" 2 "251-500" 3 ">500" label values cd4_cat cd4_cat label define hungry 1 "often" 2 "sometimes" 3 "seldom" 4 "never" label values hungry hungry label define lastalc 0 "never" 1 ">1 year ago" 2 "within last year" label values lastalc lastalc label define religion 1 "Protestant" 2 "Catholic" 3 "Saved" 4 "Moslem" 8 "Other" label values religion religion /* You can use these commands to get information about the values labels for a variable:*/ codebook religion label list religion /* Note the use of the indicator female for sex, coded 0 = no and 1 = yes. Indicators coded like this have several advantages: * the mean value of the indicator gives the proportion of women in the dataset * this is the coding required for binary outcomes in logistic models * unlike a variable for sex, coded 1 = male 2 = female, the indicator does not require us to remember the coding and value labels Now take another look at what we have. The value labels should show up in the browser. */ des sum browse /* Note that CD4 count is present only for the 30% of participants who are HIV+. ******************************************************************************** Creating and modifying new variables To create a new variable, use the generate command; to modify it, use the replace command, and to drop a variable, use the drop command. In making variables, you will often need to use so-called logical operators: equal: == not equal: ~= or ~= and: & or: | greater: > or >= less: < or <= Missing numeric values are denoted . (for missing character strings use ""). Note that missing numeric values are stored as a very large number, (always the largest in the dataset). As a result, if you are missing age values for 2 persons in the data set, and you create a new variable that includes all participants with age>30, you will include the 2 with missing data. Example: AUDIT-C is a scale measuring alcohol consumption, ranging from 0-12. For those not currently drinking (on this survey it is the past 3 months) the score is 0, for those drinking the score is 1-12. To make a new variable that represents any alcohol consumption we can use this following code: */ gen alc_any=0 if auditc==0 replace alc_any=1 if auditc>0 & auditc !=. label var alc_any "Any alcohol use" label values alc_any yesno // value label yesno is already part of the dataset * check to make sure this worked tab auditc alc_any, missing /* The second condition ensures that missings are kept as missing. A better way to do this is to use the recode command, which will handle the missing values correctly: */ drop alc_any recode auditc (0=0 "no") (2/max=1 "yes"), gen(alc_any) label var alc_any "Any alcohol use" * check to make sure this worked tab auditc alc_any, missing /* We can also use recode and replace to make a more complex indicator for unhealthy alcohol consumption, using different cutoff for men (>=4) and women (>=3): */ recode auditc (min/2=0 "no") (3/max=1 "yes"), gen(alc_unh) replace alc_unh = 0 if female==0 & auditc==3 label var alc_unh "Unhealthy alcohol use" * check to make sure this worked bysort female: tab auditc alc_unh, m /* Note the abbreviated name for the unhealthy alcohol use indicator. Stata does allow long, explicit variable names, with some limits, but they can be tedious to use in programming; better to encode the extra information in the variable label. Note also that all the alcohol variables we have made start with the same stem alc_, also with a view towards making them easier to remember. Finally, we can use recode to categorize the continuous variable, age. In doing this, we can use the ordering of the elements of the command to control whether the values at each cutpoint are put into the lower or upper category. Specifically, values equal to a cutpoint get put into the class defined by the element where the cutpoint is first encountered. We will first code age categories <=20, >20 to <=30, >30 to <=40, >40 to <=50, and >50. */ qui recode age (min/20=1 "<=20") (20/30=2 ">20 to 30") (30/40=3 ">30 to 40") /// (40/50=4 ">40 to 50") (50/max=5 ">50"), gen(age_cat) label var age_cat "Age" * check this tabstat age, by(age_cat) s(min max) /* To include the cutpoints in the upper category, we can reverse the order of the elements of the command: */ drop age_cat qui recode age (50/max=5 "50+") (40/50=4 "40 to <50") (30/40=3 "30 to <40") /// (20/30=2 "20 to <30") (min/20=1 "<20") , gen(age_cat) label var age_cat "Age" * check this tabstat age, by(age_cat) s(min max) /* Put the variables in alphabetical order, compress anything stored wastefully, get a final look at the dataet, and save the modified dataset. */ compress aorder /* Finally, save the revised dataset. It's good practice never to write over the input dataset. */ save vct_revised, replace /* Now we will use the revised dataset to do more exploration. */ use vct_revised.dta, clear /* Briefly remind yourself of what is in the dataset. Check how many observations and variables the dataset includes. */ describe /* Now get brief summary statistics for all the variables in the dataset. The full command we'll use is summarize, but we can shorten to save typing. */ sum /* Note how many people had CD4 counts. Why is this less than the number of observations? To get more information about continuous variables, add the detail option to the summarize command. */ sum age cd4, detail ******************************************************************************** /* Exploratory analysis using descriptive statistics The tabstat command gives easy summaries. Here we use the stat() option to obtain the number of observations, mean, SD, minimum, 25th percentile, median, 75th percentile, and maximum. The format option controls rounding of the statistics. */ tabstat age cd4, stat(n mean sd min p25 p50 p75 max) format(%8.3g) /* Adding the option col(statistics) transposes the results. Also, we can use the line break marker /// to tell stata that command continues after a hard return. This is handy if you want to use a relatively narrow do-file editor window. (Note that you can't use /// in the command window.) */ tabstat age cd4, stat(n mean sd min p25 p50 p75 max) /// col(statistics) format(%8.3g) * Add the by() option to see the statistics stratified by HIV status. tabstat age cd4, by(hiv) stat(n mean sd min p25 p50 p75 max) /// col(statistics) format(%8.3g) /* The tabstat command can only be used to stratify by one by variable at a time; tostratify by >1 variables, use the table command. */ table hiv female hungry, by(religion) c(mean age) format(%8.1f) * To see how categorical variables relate to each other, use the tab command. tab age_cat hiv, row tab female hiv, row tab hungry hiv, row tab alc_any hiv, row /* Taking advantage of the 0-1 coding of the HIV indicator, we can use the table command get HIV prevalence by several other variables at once; here we do just age and sex. */ table age_cat female, c(mean hiv) format(%8.2f) ******************************************************************************** /* We can also use the Stata means command to calculate confidence intervals for the means of continuous variables, both overall and within strata. */ mean age, cformat(%8.2f) mean age, over(female hiv) cformat(%8.2f) * To get, say, 90% confidence intervals, we can use the level option mean age, over(female hiv) level(90) cformat(%8.2f) /* The means command could also be used to get means and confidence intervals for binary variables coded 0-1, but a better approach is the ci prop command, with the exact option. This avoids confidence limits less than 0 or more than 1. */ ci prop hiv, exact ci prop hiv, exact level(90) /* The ci command does not allow stratification, but we can get around that using the bysort prefix: */ bysort female: ci prop hiv, exact * The ci command also works for continous and Poisson-distributed count varibles ci means age ci means n_adults, poisson ******************************************************************************** /* Exploratory data analysis using graphical means. First make a histogram of the age distribution. We can use the name option so that multiple graphs can be shown in the viewer. The replace option is needed to prevent Stata from stopping if the graph already exists -- again, this is useful if you are rerunning a do-file to debug it. */ histogram age, fcolor(blue) name(hist_age_00, replace) /* Now add titles to the histogram. */ histogram age, fcolor(blue) title(Histogram of age) /// xtitle(Age) name(hist_age_01, replace) /* By default, the Y-axis shows the "density", which sums to one for all bars. Alternatively, we could get a histogram of frequencies. */ histogram age, fcolor(blue) frequency /// title(Histogram of age) xtitle(Age) name(hist_age_02, replace) /* Now combine the two histograms into a single plot. Note that he histograms are identical in shape, and both show that age is somewhat skewed to the right. */ graph combine hist_age_01 hist_age_02, name(hist_age_combined, replace) /* To save this as a PDF file, we can use the following command. Other options include saving as png, wmf, tif, eps (encapsulated postscript), and svg (scalable vector graphics), which some journals prefer. See help graph export for more information. */ graph export hist_age_combined.pdf, replace /* Now we will change the number of bars in the histogram, leaving Stata to decide on their widths. We will do this using "loops." Specifically, the program will loop over the 4 values of nbins (i.e., # of bins) in the first line, pasting those values as so-called local variables `nbins' into the histogram command in several places. Learning to use loops can make programming a lot less tedious. You will need to run the four lines of code defining the loop together. */ foreach nbins in 5 10 25 50 { histogram age, fcolor(blue) percent bin(`nbins') /// title("`nbins' bins") xtitle(Age) name(hanb`nbins', replace) } graph combine hanb5 hanb10 hanb25 hanb50, name(hist_age_bins, replace) /* As an alternative, we can vary the width of the bins, again using a loop. */ foreach bw in 1 2 5 10 { histogram age, fcolor(blue) percent width(`bw') /// title(bin width `bw') xtitle(Age) name(habw`bw', replace) } graph combine habw1 habw2 habw5 habw10, name(hist_age_widths, replace) ******************************************************************************** /* Box plots are also useful for assessing the distribution of continuous variables. The box runs from the 25th to the 75th percentile, and includes a line at the median. The so-called whiskers extend 1.5 IQRs from the bottom and top of the box, or to the most extreme value, whichever is closer to the box. Any outliers more than 1.5 IQRs from the box are plotted separately. */ graph box age, name(box_age, replace) /* The boxplot also shows that age is somewhat right skewed, as we saw in the histograms. We can stratify the boxplot by up to three other factors. Note that the value labels we added in lab 1 to sex and hiv make this plot easier to understand */ graph box age, over(sex) over(hiv) name(box_age_stratified, replace) /* Finally, Q-Q (quantile-quantile) plots are a good way to assess Normality. We will use a loop to look at them within groups, using some new commands to control the loops and title the plots: */ qui levelsof sex, local(sl) // extracts distinct values of sex foreach i in `sl' { local sex : label(sex) `i' // gets value label for level i of sex forvalues j = 0/1 { local stat : label(hiv) `j' // gets value label for level j of hiv qnorm age if sex==`i' & hiv==`j', /// title("`stat' `sex'") name(qq_age`i'`j', replace) } } graph combine qq_age10 qq_age11 qq_age20 qq_age21, name(qq_age, replace) /* With the exception of HIV positive men, these plots all show the same right skewness we observed in the box plots and histograms. */ ******************************************************************************** /* Now we will use scatter plots to look at the association of a continuous outcome with a continuous predictor. */ scatter n_support age, name(scatter_age, replace) /* Now we add a Lowess smooth through the data. Essentially this is a regression of n_support on age, without making the assumption that the regression line is linear. The bwidth option controls the smoothness of the line. The default is 0.8, which means that Stata uses the nearest 80% of the data to determine the vertical location of the smooth at each value of the independent variable. The lineopts() option allows us to control the thickness of the line, and could also be used to reset its color and whether the line is solid, dashed, etc. We will again using loops to explore the bandwidth. The looping variable is a percentage, because the plot name can't include decimals. This means that we have to use some arcane code to convert `pct' to a proportion for use in the bw() option in the lowess command. Note that the first line of the loop uses a somewhat different syntax than the previous loops. */ foreach pct of numlist 80(-20)20 10 1 { lowess n_support age, /// bwidth(`=`pct'/100') msize(tiny) lineopts(lwidth(medthick)) title("") /// name(bw`pct', replace) } graph combine bw80 bw60 bw40 bw20 bw10 bw1, name(lowess_combined, replace) /* Reducing the bandwith imposes less smoothness, but clearly the smallest bandwidth is going too far! We can also focus on a certain age range by adding an if statement. */ lowess n_support age if age<60, /// msize(tiny) lineopts(lwidth(medthick)) title("") /// name(lowess_restricted, replace) /* Finally, we can use tabstat and box plots to look at the conditional distribution of a continuous outcome when the independent variable is categorical. */ tabstat n_support, by(hungry) stat(n mean sd min p25 p50 p75 max) format(%8.3g) graph box n_support, over(hungry) name(ins_alc, replace) /* Now close the log file.*/ log close