* Feb 2, 2005, LAB 3
* Read in the three cows data sets and create a STATA data set
***First open a .log file to make sure you keep all you're doing for later use. The .do file has only ***commands, while the .log file also keeps STATA output, with the exception of graphs.
***Make sure you change the path and replace it with your own 
log using "C:\My Documents\655LDA\DataLDA\cows-readin.log" , replace
set memory 40m
set matsize 100
* read in three data sets, starting with last. The file is space separated so wee use the infile command

infile prot1-prot19 using "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.lupins.data"

gen str6 diet="lupins"
compress
***Now save the STATA file created in a location you want 
save  "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.dta" , replace
clear

infile prot1-prot19 using "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.mixed.data"
gen str6 diet="mixed"
compress
append using  "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.dta"
save  "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.dta", replace
clear

infile prot1-prot19 using "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.barley.data"
gen str6 diet="barley"
compress
append using  "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.dta" 

order diet
compress

* Set missing values and save in the wide format 

for X in num 1/19 : replace protX=. if (protX==0.00)
save "C:\Documents and Settings\Sorina\Desktop\LDA\Lab3\cows.dta" , replace


* now make in long format and save
gen id = _n  
*this says that the id is the observation number (n)
order id
compress
reshape long prot, i(id) j(week)

order id diet week 
**label the variables so you know what they stand for. This is completely optional. 
label var id   "Cow ID"
label var diet "Diet"
label var week "Week since calving"
label var prot "Protein content"

*Spaghetti plots

sort diet id week
***STATA 8 version
twoway line prot week, by(diet) c(L) s(.) 
**STATA 7 version of the above commans is     gr prot week, by(diet) c(l) s(i)

*Kernel and smoothing
 ksm prot week, lowess gen(cowsmth)
drop cowsmth ***Drops a variable that I don't want to use later 

lowess prot week if diet=="barley", gen(sm1) bw(.6) nograph
lowess prot week if diet=="mixed", gen(sm2) bw(.6) nograph
lowess prot week if diet=="lupins", gen(sm3) bw(.6) nograph
label var sm1 "barley"
label var sm2 "mixed"
label var sm3 "lupins"

sort week
twoway line sm1 sm2 sm3 week, c(lll) s(x) xlab ylab l1("Protein content") 

drop sm1-sm3
****Convert an ordinary dataset into a longitudinal dataset (time series)
tsset id week
**Alternatively, one could use the iis and tis commands
*iis id
*tis week

**Exploratory Data Analysis
describe
tsreport , report panel
xtdes
xtsum prot

 xtgraph prot, group(diet) ti("Mean protein content of milk") bar(no)

*Estimating variance within subjects and between subjects 
xtsumcorr prot
anova prot id
*Make a Variogram without removing covariate effects
variogram protrs , discrete incl(100)

*Make a Variogram removing covariate effects

* Compute residuals, removing effects of diet and time
sort diet week
by diet week : egen protmn=mean(prot)
gen protrs = prot - protmn
sort id week
xtsumcorr protrs

* Compute discrete-lag variogram
variogram protrs , discrete incl(100)
autocor protrs week id

save   "C:\My Documents\655LDA\DataLDA\cows.dta" , replace
clear
log close