/* Biostat 200 Lab 1b capture log close log using "lab 1b", replace /* Read in the data posted on the class web site. Again, these data come from a study of voluntary counseling and testing (VCT) for HIV at Mulago Hospital in Kampala, Uganda. This is the version of the dataset we worked on last week. */ 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) ******************************************************************************** /* 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