* Stata code for 2010 ABCT workshop on HLM for Longitudinal Data Analysis
clear
set more off
*read in data
insheet using "Sex and Stress Data UPDATED.txt", tab
*take a look at data structure
describe, detail
*first 10 rows of data
list in 1/10
*create a numeric gender variable (men=1, women=0)
gen gender2=gender=="Male"
*summary of variables
summarize week-sex stress-gender2, detail
*sort by id and week
sort id week
*Is data balanced by week and id?
tabulate week
tabulate id
*What does outcome look like?
histogram sex, freq discrete xtitle("") ytitle("Frequency") title("Histogram of Number of Sexual Encounters (per week)") text(250 20 "M: 2.6") text(230 20 "SD: 2.9")
*Growth curves for 20 individuals
twoway line sex week if (id>=75 & id<=94), by(id, note("") title("Number of Sexual Encounters" "Across Weeks for 20 Subjects")) xtitle("Week") ytitle("Sex") xlabel(0(2)12) ylabel(0(3)15)
*Spaghetti Plot
spagplot sex week, id(id) nofit by(gender)
*Fit HLMs
xtmixed sex || id:, var
*grand mean center stress
quietly summarize stress
generate stress_c=stress-r(mean)
*stress main-effects model
xtmixed sex stress_c || id:,var
*add gender main effect
xtmixed sex stress_c gender2 || id:, var
*add interaction
xtmixed sex c.stress_c##i.gender2 || id:, var
estimates store mod1
estat ic
*do we need a random effect for stress?
xtmixed sex c.stress_c##i.gender2 || id:stress_c, var cov(un)
estimates store mod2
estat ic
lrtest mod2 mod1
*is there any residual "week" effect?
xtmixed sex c.stress_c##i.gender2 week || id:stress_c, var cov(un)
estat ic
*create residuals,random effects, and fitted values
predict residuals, residuals
predict reffect*, reffects
predict fitted, fitted
*histograms and q-q plots of residuals and random effects
histogram(residuals)
qnorm(residuals)
histogram(reffect1)
qnorm(reffect1)
histogram(reffect2)
qnorm(reffect2)
*scatterplot of residuals against fitted values
twoway scatter resid fitted, jitter(2)
*calculating means of stress and sex within person
*can use these to calculate centered values
by id, sort: egen stressmn=mean(stress)
by id, sort: egen sexmn=mean(sex)
*couple data
clear
insheet using "Atkins 2005 wide w missing.txt"
destring das, replace
*unconditional model
xtmixed das || id: || sex: , var
estimates store model1
*add time and therapy
xtmixed das therapy time || id: || sex: , var
estimates store model2
*add time random-effect at couple-level
xtmixed das therapy time || id:time || sex:, var
estimates store model3
*add time random-effect at individual-level
xtmixed das therapy time || id:time || sex:time, var
发帖
雷达卡
京公网安备 11010802022788号
