* Get the starting values for PROC NLMIXED; 
PROC MIXED DATA = health_behav2 METHOD = ml COVTEST;
CLASS id;
MODEL lga = sq/ SOLUTION;
RANDOM intercept/TYPE=un SUBJECT =id;
REPEATED/LOCAL = EXP(sq);
ODS OUTPUT SolutionF = beta;
ODS OUTPUT CovParms = random;
RUN;

* MELS Model in which all subjects are treated equally;
* Results used as starting values for leave-one-out MELS models; 
PROC NLMIXED DATA = health_behav2 TECH = quanew GCONV=1e-12 METHOD = gauss 
COV ECOV;
PARMS beta0 = 1.8795 beta_sq = 0.1038
tau0 =  -0.2172857 tau_sq = -0.00898
alpha0 = -1.129793
StdDev_scale = 0.05
corr = 0.05;

varBS = EXP(alpha0);
nu = SQRT(varBS) * theta1;

omega = StdDev_scale * theta2;

mu = beta0  + beta_sq * sq + nu;
varWS = EXP(tau0 + tau_sq * sq + omega);
MODEL lga ~ NORMAL(mu, varWS);
RANDOM theta1 theta2 ~ NORMAL([0,0], [1, corr, 1]) SUBJECT = id;

* Get the estimates for variances and covariance of random effects;
ESTIMATE "varBS" EXP(alpha0);
ESTIMATE "varWS" StdDev_scale**2; 
ESTIMATE "COV" corr * StdDev_scale * SQRT(EXP(alpha0)); 

ODS OUTPUT ParameterEstimates = healthbehav_all_parm;
ODS OUTPUT FitStatistics= healthbehav_all_fit;
ODS OUTPUT CovMatParmEst= healthbehav_all_cov;
ODS OUTPUT AdditionalEstimates= healthbehav_all_est;
ODS OUTPUT CovMatAddEst= healthbehav_all_est_cov;
    
RUN;


* Get all the distinct IDs in the data set; 
PROC SQL NOPRINT;
  SELECT DISTINCT id INTO :unique_id separated by ' ' FROM health_behav2;
QUIT;

* Re-estimate the model with each subject SEPARATED from the estimation of; 
* both random location effects and random scale effects;
%MACRO INFLUENTIAL;
/*store:
1. sub_parm: fixed effect estimates
2. sub_fit: model fit statistics (-2 log likelihood, AIC, AICC, BIC)
3. sub_cov: covariance matrix of fixed effect estimates
4. sub_est: estimates of variances and covariance of random effects
5. sub_est_cov: variance-covariance matrix of the variances and 
covariance of random effects*/
%DO i = 1 %TO %SYSFUNC(COUNTW(&unique_id));
%LET this_id = %SCAN(&unique_id,&i,%STR( ));
	ODS OUTPUT ParameterEstimates = outlier.sub_parm_new;
	ODS OUTPUT FitStatistics = outlier.sub_fit_new;
	ODS OUTPUT CovMatParmEst = outlier.sub_cov_new;
	ODS OUTPUT AdditionalEstimates= outlier.sub_est_new;
	ODS OUTPUT CovMatAddEst= outlier.sub_est_cov_new;

	PROC NLMIXED DATA = r2_sim.health_behav2 TECH = quanew GCONV=1e-12 
 METHOD = gauss COV ECOV;
	PARMS beta0 = 1.8795 beta_sq = 0.1038
	tau0 =  -0.2172857 tau_sq = -0.00898
	alpha0 = -1.129793
	StdDev_scale = 0.0001
	c = 0.0001
	d = 0.0001
	corr = 0.0001;

	varBS = EXP(alpha0);
	nu = SQRT(varBS) * theta1;
	
	omega = StdDev_scale * theta2;
	
	mu = beta0 * (id ~= &this_id) + beta_sq * sq + nu * (id ~= &this_id) + 
 c * (id = &this_id);
	varWS = EXP(tau0 * (id ~= &this_id) + tau_sq * sq + 
 omega * (id ~= &this_id) + d * (id = &this_id));
	MODEL lga ~ NORMAL(mu, varWS);
	RANDOM theta1 theta2 ~ NORMAL([0,0], [1, corr, 1]) SUBJECT = id;
	ESTIMATE "varBS" EXP(alpha0);
	ESTIMATE "varWS" StdDev_scale**2; 
	ESTIMATE "COV" corr * StdDev_scale * SQRT(EXP(alpha0)); 
	
	RUN; 
	
	* Select variables and add subject ID;
	DATA outlier.sub_parm_new;
       SET outlier.sub_parm_new;
       id=&this_id;
	RUN;
	DATA outlier.sub_est_new;
       SET outlier.sub_est_new;
       id=&this_id;
	RUN;
	DATA outlier.sub_fit_new;
       SET outlier.sub_fit_new;
       id=&this_id;
       KEEP Descr Value ID;
	RUN;
	DATA outlier.sub_cov_new;
       SET outlier.sub_cov_new;
       id=&this_id;
       DROP Row;
	RUN;
	DATA outlier.sub_est_cov_new;
       SET outlier.sub_est_cov_new;
       id=&this_id;
       DROP Row;
	RUN;
	
PROC APPEND BASE=outlier.sub_parm NEW=outlier.sub_parm_new;
RUN;
PROC APPEND BASE=outlier.sub_est NEW=outlier.sub_est_new;
RUN;
PROC APPEND BASE=outlier.sub_fit NEW=outlier.sub_fit_new;
RUN;
PROC APPEND BASE=outlier.sub_cov NEW=outlier.sub_cov_new;
RUN;
PROC APPEND BASE=outlier.sub_est_cov NEW=outlier.sub_est_cov_new;
RUN;

%END; 
%MEND INFLUENTIAL;

%INFLUENTIAL;