/* Generated Code (IMPORT) */
/* Source File: SAS_Data_A_series.xlsx */
/* Source Path: /home/u44778551/Jan Biel */
/* Code generated on: 1/17/20, 6:41 AM */

%web_drop_table(WORK.IMPORT);


FILENAME REFFILE '/home/u56336265/Original data/20230709 A-C 20-40-80 CS20 US5X100 unpaired.xlsx';

PROC IMPORT DATAFILE=REFFILE
	DBMS=XLSX
	OUT=WORK.IMPORT;
	GETNAMES=YES;
RUN;

PROC CONTENTS DATA=WORK.IMPORT; RUN;

Title1 '20230709 A-C 20-40-80 CS20 US5X100 unpaired';


%web_open_table(WORK.IMPORT);

/*
PROC PRINT;
RUN;
*/

data Eyes;
set Work.Import;

if Intval =< 0 then delete;
Freq = 1000/Intval; /* Frequency in Hz */

AdjMidTime = MidTime-InitStim;

   proc sort data=eyes;
   by Series Contrast Width Speed Rep;
   run;

/*
   proc print;
   Title '20-40-80 A  and C series';
   run;
*/  
   

/* Check data */
proc means mean min max data=eyes noprint;
var Freq;
output out= averages mean= MeanFrequency min = MinFrequency max = MaxFrequency;
by  Series Contrast Width Speed;
run;

/*
proc print data=averages;
var Series Contrast Width Speed MeanFrequency MinFrequency MaxFrequency ;
run;
*/

/* replace original values with residuals */
proc glm data= eyes noprint;
/* Class Speed Stimulus Width Eye; */
model Freq =  /solution;
output out= residuals r= res;
by Series Contrast Width Speed;
run;

/*
proc print data= residuals;
run;
*/


proc sort data=Residuals;
by Series Contrast Width Speed AdjMidTime Res;
run;

/*
proc print data=residuals;
var Series Contrast Width Speed AdjMidTime Res;
run;
*/

title2 'Extreme Observations';
ods select ExtremeObs;
proc univariate data=Residuals;
   var Res;
   id AdjMidTime;
by Series Contrast Width Speed;
run;

  /* Model average response of treatments*/
  data Adjusted;
  set Residuals;

  proc sort data=adjusted;
  title2 'Relationship between time and neural activity';
  by Contrast Width Speed AdjMidTime;
  run;
 
  ods rtf file = '/home/u56336265/Output data/20230709 A-C 20-40-80 CS20 US5X100 unpaired.rtf';
  
  proc nlin data=adjusted noitprint alpha = 0.05 ;
  Title2 '(full model)';
  parms 
  Level = 0
  Reaction = 5
  difLevel = 0
  difPeak = 0
  difa = 0
  difb = 0
  a = 500
  b = 500;
  bounds Reaction >= 0;

  Time1 =  -20000;
  Time2 = 30000;
  
  if Series = 'A' then Level1 = Level; else Level1 = Level + DifLevel;
  if Series = 'A' then Peak = Reaction; else Peak = Reaction + DifPeak;
  if Series = 'A' then a1 = a; else a1 = a + Difa;
  if Series = 'A' then b1 = b; else b1 = b + Difb;
    
  if (AdjMidTime > Time1) and (AdjMidTime < Time2) then F = exp(-a1*log(a1/(a1+b1))- b1*log(b1/(a1+b1)) + log(Peak) + a1*log((AdjMidTime-Time1)/(Time2-Time1)) + b1*log((Time2-AdjMidTime)/(Time2-Time1)));
  else F = 0; 
    
  model Res = Level1 + F; 
  output out= nlinCurves pred= PredFrequency l95m= lower u95m= upper student = ttest;
  by Contrast Width Speed;
  run;
  
  proc nlin data=adjusted noitprint alpha = 0.05 ;
  Title2 '(reduced model)';
  parms 
  Level = 0
  Reaction = 5
  difLevel = 0
  difPeak = 0
  a = 500
  b = 500;
  bounds Reaction >= 0;

  Time1 =  -30000;
  Time2 = 50000;
  
  if Series = 'A' then Level1 = Level; else Level1 = Level + DifLevel;
  if Series = 'A' then Peak = Reaction; else Peak = Reaction + DifPeak;
      
  if (AdjMidTime > Time1) and (AdjMidTime < Time2) then F = exp(-a*log(a/(a+b))- b*log(b/(a+b)) + log(Peak) + a*log((AdjMidTime-Time1)/(Time2-Time1)) + b*log((Time2-AdjMidTime)/(Time2-Time1)));
  else F = 0; 
  
  
  model Res = Level1 + F; 
  output out= nlinCurves pred= PredFrequency l95m= lower u95m= upper student = ttest;
  by Contrast Width Speed;
  run;

ods rtf close; 
 
 
  
  