/* Generated Code (IMPORT) */
/* Source File: Original data.xlsx */
/* Source Path: /home/u56336265/Associative learning SAS scripts */
/* Code generated on: 1/17/20, 6:41 AM */

%web_drop_table(WORK.IMPORT);


FILENAME REFFILE '/home/u56336265/Original data/20230709 A5 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 A5 20-40-80 CS20 US5X100 unpaired R2';


%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 Contrast Width Speed Rep;
run;

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

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

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

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

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

  proc sort data=adjusted;
  by Contrast Width Speed AdjMidTime;
  run;

   
  ods rtf file = '/home/u56336265/Output data/20230709 A5 20-40-80 CS20 US5X100 unpaired R2 (adjusted).rtf';
    
  proc nlin data=adjusted noitprint ;
  title2 'Relationship time and neural activity';
  parms 
  Level = 0
  Reaction = 5
  a = 500
  b = 500;
  bounds Reaction >= 0;

  Time1 =  -30000;
  Time2 = 50000;
  if (AdjMidTime > Time1) and (AdjMidTime < Time2) then F = exp(-a*log(a/(a+b))- b*log(b/(a+b)) + log(Reaction) + a*log((AdjMidTime-Time1)/(Time2-Time1)) + b*log((Time2-AdjMidTime)/(Time2-Time1)));
  else F = 0; 
    
  model Res = Level + F; 
  output out= nlinCurves pred= PredFrequency l95m= lower u95m= upper;
  by Contrast Width Speed;
  run;
  
  /*proc print data = nlinCurves;
  run;*/
  

  proc sort data= nlinCurves;
  by Contrast Width Speed;
  run;


 
  proc gplot data=nlinCurves;
  title2 'Adjusted nlinCurves';
  plot Res*AdjMidTime PredFrequency*AdjMidTime lower*AdjMidTime Upper*AdjMidTime/ overlay haxis= -30000 to 50000 by 10000;
  axis1 length=70000;
  symbol1 interpol =  none value= plus cv=blue;
  symbol2 interpol=join ci= red width= 3;
  symbol3 interpol=join ci= black width= 1;
  symbol4 interpol=join ci= black width= 1;
  by Contrast Width Speed; 
  run;
  
   proc gplot data=nlinCurves;
  title2 'Adjusted nlinCurves';
  plot  Res*AdjMidTime PredFrequency*AdjMidTime lower*AdjMidTime Upper*AdjMidTime/ overlay haxis= -20000 to 20000 by 10000;
  axis1 length=70000;
  symbol1 interpol =  none value= plus cv=blue;
  symbol2 interpol=join ci= red width= 3;
  symbol3 interpol=join ci= black width= 1;
  symbol4 interpol=join ci= black width= 1;
  by Contrast Width Speed; 
  run;
  

  
/* Compute R2 */
data explained;
set nlinCurves;

/* only positive values of predicted frequencies are used in the calculations of R2 */
if PredFrequency < 0 then delete;

proc reg data= Explained;
title 'Model fit';
model Res = PredFrequency;
by Contrast Width Speed;
run;

ods rtf close;
   
/* export data as an Excel file */ 
PROC EXPORT 
  data = nlinCurves
  dbms =xlsx
  outfile = "/home/u56336265/Output nlin/20230709 A5 20-40-80 CS20 US5X100 unpaired R2 nlin curves (adjusted)"
  REPLACE;
RUN;
  
  