# First the legal disclaimer
# Program name: SemProgramNing.txt
# Copyright Yale Pepper Center (2008)  author: Yuming Ning 

   
#   This program is free software: you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation, either version 3 of the License, or
#   (at your option) any later version.
   
#   This program is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
   
#   You should have received a copy of the GNU General Public License
#   along with this program.  If not, see <http://www.gnu.org/licenses/>.


#The following Mplus code is used to analyze a structural equation model (SEM) for 
#the relationship among the dischargement from skilled-nursing facilities, the situation for
#activities of daily living (ADL), the geriatric conditions, and the comorbidities of chronic
#diseases in a SNFPAC study. The SNFPAC study is a prospective cohort study of Medicare
#beneficiaries admitted for post-acute rehab to a large skilled-nursing facility in
#New Haven, CT. 
#The hypothesized pathways of the SEM is as follows:
#1. number of geriatric conditions-->latent ADL effect-->discharge back home;
#2. latent ADL score contains initial ADL & ADL change during study time;
#3. number of chronic comorbidities-->discharge back home;
#4. discharge back home is controlled by gender, medicaid availability, and orthopedic treat;
 

############################# CODE BEGIN ############################### 



/*To Remove an observed value from a dataset to simulate a missingness, and then 
  to impute the pseudo 'missing value' with different multiple imputation methods
  and compare the imputation results*/ 


*Bulid the original working dataset with 7 waves of observations; 
data chs1.origin;
merge chs.baseline chs.followup_2 chs.followup_3 chs.followup_4 chs.followup_5 chs.followup_6 chs.followup_7;
by idno;
run;


*Record real values of observed data that will be removed later;
data chs1.newadl; 
set chs1.origin; 
array newadl [7] newadl1-newadl7;
do j=1 to 4;
 if newadl[j]>=0 & newadl[j+1]<0 & newadl[j+2]>=0 & newadl[j+3]>=0 then do;
  value=newadl[j+2];
  time=j+2;
  keep idno time value;
  output;
 end;
end;
run;


*Remove some observed values of the newadl variable from original data to simulate missingness;
data chs1.simumiss;
set chs1.origin;
array newadl [7] newadl1-newadl7;

do j=1 to 4;
 if newadl[j]>=0 & newadl[j+1]<0 & newadl[j+2]>=0 & newadl[j+3]>=0 then do;
  newadl[j+2]=.;
  t_adl=j+2;
 end;
end;

drop j;
run;


*Generate 7 datasets with pseudo missing newadl values for 7 waves of observations;
*which will be used for the sequential multiple imputation method;
data chs1.simumiss1;
set chs1.simumiss;
keep idno perstat agebl female nonwhite 
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1;
run;
data chs1.simumiss2;
set chs1.simumiss;
keep idno perstat agebl female nonwhite 
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1
     newadl2 newgenhel2 newgrip2 newstand2 newrise2 newtime2 
     new_wt2 newstep2 newmile2 newmmse2 newdig2 socsup2 newdep2
	 newhear2 newvision2;
run;
data chs1.simumiss3;
set chs1.simumiss;
keep idno perstat agebl female nonwhite
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1
     newadl2 newgenhel2 newgrip2 newstand2 newrise2 newtime2 
     new_wt2 newstep2 newmile2 newmmse2 newdig2 socsup2 newdep2
	 newhear2 newvision2
     newadl3 newgenhel3 newgrip3 newstand3 newrise3 newtime3 
     new_wt3 newstep3 newmile3 newmmse3 newdig3 socsup3 newdep3
	 newhear3 newvision3;
run;
data chs1.simumiss4;
set chs1.simumiss;
keep idno perstat agebl female nonwhite
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1
     newadl2 newgenhel2 newgrip2 newstand2 newrise2 newtime2 
     new_wt2 newstep2 newmile2 newmmse2 newdig2 socsup2 newdep2
	 newhear2 newvision2
     newadl3 newgenhel3 newgrip3 newstand3 newrise3 newtime3 
     new_wt3 newstep3 newmile3 newmmse3 newdig3 socsup3 newdep3
	 newhear3 newvision3
     newadl4 newgenhel4 newgrip4 newstand4 newrise4 newtime4 
     new_wt4 newstep4 newmile4 newmmse4 newdig4 socsup4 newdep4
	 newhear4 newvision4;
run;
data chs1.simumiss5;
set chs1.simumiss;
keep idno perstat agebl female nonwhite
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1
     newadl2 newgenhel2 newgrip2 newstand2 newrise2 newtime2 
     new_wt2 newstep2 newmile2 newmmse2 newdig2 socsup2 newdep2
	 newhear2 newvision2
     newadl3 newgenhel3 newgrip3 newstand3 newrise3 newtime3 
     new_wt3 newstep3 newmile3 newmmse3 newdig3 socsup3 newdep3
	 newhear3 newvision3
     newadl4 newgenhel4 newgrip4 newstand4 newrise4 newtime4 
     new_wt4 newstep4 newmile4 newmmse4 newdig4 socsup4 newdep4
	 newhear4 newvision4
     newadl5 newgenhel5 newgrip5 newstand5 newrise5 newtime5 
     new_wt5 newstep5 newmile5 newmmse5 newdig5 socsup5 newdep5
	 newhear5 newvision5;
run;
data chs1.simumiss6;
set chs1.simumiss;
keep idno perstat agebl female nonwhite
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1
     newadl2 newgenhel2 newgrip2 newstand2 newrise2 newtime2 
     new_wt2 newstep2 newmile2 newmmse2 newdig2 socsup2 newdep2
	 newhear2 newvision2
     newadl3 newgenhel3 newgrip3 newstand3 newrise3 newtime3 
     new_wt3 newstep3 newmile3 newmmse3 newdig3 socsup3 newdep3
	 newhear3 newvision3
     newadl4 newgenhel4 newgrip4 newstand4 newrise4 newtime4 
     new_wt4 newstep4 newmile4 newmmse4 newdig4 socsup4 newdep4
	 newhear4 newvision4
     newadl5 newgenhel5 newgrip5 newstand5 newrise5 newtime5 
     new_wt5 newstep5 newmile5 newmmse5 newdig5 socsup5 newdep5
	 newhear5 newvision5
     newadl6 newgenhel6 newgrip6 newstand6 newrise6 newtime6 
     new_wt6 newstep6 newmile6 newmmse6 newdig6 socsup6 newdep6
	 newhear6 newvision6;
run;
data chs1.simumiss7;
set chs1.simumiss;
keep idno perstat agebl female nonwhite
     newadl1 newgenhel1 newgrip1 newstand1 newrise1 newtime1 
     new_wt1 newstep1 newmile1 newmmse1 newdig1 socsup1 newdep1
	 newhear1 newvision1
     newadl2 newgenhel2 newgrip2 newstand2 newrise2 newtime2 
     new_wt2 newstep2 newmile2 newmmse2 newdig2 socsup2 newdep2
	 newhear2 newvision2
     newadl3 newgenhel3 newgrip3 newstand3 newrise3 newtime3 
     new_wt3 newstep3 newmile3 newmmse3 newdig3 socsup3 newdep3
	 newhear3 newvision3
     newadl4 newgenhel4 newgrip4 newstand4 newrise4 newtime4 
     new_wt4 newstep4 newmile4 newmmse4 newdig4 socsup4 newdep4
	 newhear4 newvision4
     newadl5 newgenhel5 newgrip5 newstand5 newrise5 newtime5 
     new_wt5 newstep5 newmile5 newmmse5 newdig5 socsup5 newdep5
	 newhear5 newvision5
     newadl6 newgenhel6 newgrip6 newstand6 newrise6 newtime6 
     new_wt6 newstep6 newmile6 newmmse6 newdig6 socsup6 newdep6
	 newhear6 newvision6
     newadl7 newgenhel7 newgrip7 newstand7 newrise7 newtime7 
     new_wt7 newstep7 newmile7 newmmse7 newdig7 socsup7 newdep7
	 newhear7 newvision7;
run;


*********************************************************
Now, go to use PROC MI to impute the simulated datasets
and get the imputed data back here
*********************************************************;


*Evaluate the sequential multiple imputations of newadl;
data chs1.adl_se1;
merge chs1.newadl (in=a) chs1.impse1_7i_1;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se2;
merge chs1.newadl (in=a) chs1.impse1_7i_2;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se3;
merge chs1.newadl (in=a) chs1.impse1_7i_3;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se4;
merge chs1.newadl (in=a) chs1.impse1_7i_4;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se5;
merge chs1.newadl (in=a) chs1.impse1_7i_5;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se6;
merge chs1.newadl (in=a) chs1.impse1_7i_6;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se7;
merge chs1.newadl (in=a) chs1.impse1_7i_7;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se8;
merge chs1.newadl (in=a) chs1.impse1_7i_8;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se9;
merge chs1.newadl (in=a) chs1.impse1_7i_9;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_se10;
merge chs1.newadl (in=a) chs1.impse1_7i_10;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;


*Evaluate the simultaneous multiple imputations of newadl;
data chs1.adl_si1;
merge chs1.newadl (in=a) chs1.impsi1_7_1;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si2;
merge chs1.newadl (in=a) chs1.impsi1_7_2;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si3;
merge chs1.newadl (in=a) chs1.impsi1_7_3;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si4;
merge chs1.newadl (in=a) chs1.impsi1_7_4;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si5;
merge chs1.newadl (in=a) chs1.impsi1_7_5;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si6;
merge chs1.newadl (in=a) chs1.impsi1_7_6;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si7;
merge chs1.newadl (in=a) chs1.impsi1_7_7;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si8;
merge chs1.newadl (in=a) chs1.impsi1_7_8;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si9;
merge chs1.newadl (in=a) chs1.impsi1_7_9;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;
data chs1.adl_si10;
merge chs1.newadl (in=a) chs1.impsi1_7_10;
by idno;
if a;
array newadl [7] newadl1-newadl7;
do i=1 to 7;
 if i=time then do;
  value1=newadl[i];
  time1=i;
  diff=value1-value;
  diffa=abs(value1-value);
  diff2=(value1-value)**2;
  time_diff=time1-time;
 end;
end;
keep idno value value1 diff diffa diff2 time time1 time_diff;
run;