Merge pull request #186 from drbenvincent/dev

Dev

.gitignore

@@ -6,15 +6,13 @@
 
 figs/
 
-
 demo/figs/*
 demo/output/*
 
 */temp.mat
 
 */COMBINED_DATA.txt
 
-
 demo/data/methodspaper-kirby27.txt
 
 demo/parameterEstimates_methodspaper-kirby27.txt
@@ -23,8 +21,6 @@ demo/parameterEstimates_nonHierarchical.txt
 
 demo/*.csv
 
-
-
 ddToolbox/models/hierarchicalME
 
 ddToolbox/models/hierarchicalLogK
@@ -53,13 +49,12 @@ output-2.csv
 
 temp.data.R
 
-demo/datasets/non-parametric/groupLevelData/COMBINED_DATA.txt
+demo/datasets/non_parametric/groupLevelData/COMBINED_DATA.txt
 
 demo/datasets/kirby/groupLevelData/COMBINED_DATA.txt
 
 tests/tests/grw_test/ConvergenceReport.txt
 
-
 output/*
 
 demo/datasets/test_data/groupLevelData/COMBINED_DATA.txt
@@ -113,3 +108,7 @@ ddToolbox/models/parametric_models/hyperbolic_magnitude_effect_models/hierarchic
 ddToolbox/models/parametric_models/hyperbolic_models/separateLogK
 
 demo/datasets/non_parametric/groupLevelData/COMBINED_DATA.txt
+
+ddToolbox/models/parametric_models/hyperbolic_models/hierarchicalLogK
+
+demo/temp.init.R

README.md

@@ -1,6 +1,6 @@
 # Hierarchical Bayesian estimation and hypothesis testing for delay discounting tasks
 
-Vincent, B., T. (in press) **[Hierarchical Bayesian estimation and hypothesis testing for delay discounting tasks](http://link.springer.com/article/10.3758%2Fs13428-015-0672-2)**, Behavior Research Methods. doi:10.3758/s13428-015-0672-2
+Vincent, B. T. (2016) **[Hierarchical Bayesian estimation and hypothesis testing for delay discounting tasks](http://link.springer.com/article/10.3758%2Fs13428-015-0672-2)**, Behavior Research Methods. 48(4), 1608-1620. doi:10.3758/s13428-015-0672-2
 
 [![Issue Stats](http://issuestats.com/github/drbenvincent/delay-discounting-analysis/badge/issue?style=flat-square)](http://issuestats.com/github/drbenvincent/delay-discounting-analysis)
 
@@ -14,9 +14,9 @@ Key features:
 * Parameters exported to a `.csv` file for analysis in [JASP](https://jasp-stats.org).
 * Optionally use hierarchical inference to improve participant-level estimates.
 * A variety of models are available:
-  * Exponential discounting.
-  * Hyperbolic discounting.
-  * Hyperbolic discounting + magnitude effect, where discount rates vary as a function of reward magnitude.
+  * 1-parameter discount functions: exponential, hyperbolic.
+  * 2-parameter discount functions: hyperboloid
+  * Also, hyperbolic discounting + magnitude effect, where discount rates vary as a function of reward magnitude.
 * Explicit modelling of participant errors provides more robust parameter estimates of discounting parameters.
 * Posterior predictive checks help evaluate model goodness and aid data exclusion decisions.
 * Publication quality figures.
@@ -38,7 +38,7 @@ The commands use to get your analysis up and running are quite quick and easy. A
 ## More information
 
 ### Read the research paper
-Vincent, B., T. (in press) [Hierarchical Bayesian estimation and hypothesis testing for delay discounting tasks](http://link.springer.com/article/10.3758%2Fs13428-015-0672-2), Behavior Research Methods. doi:10.3758/s13428-015-0672-2
+Vincent, B. T. (2016) [Hierarchical Bayesian estimation and hypothesis testing for delay discounting tasks](http://link.springer.com/article/10.3758%2Fs13428-015-0672-2), Behavior Research Methods. 48(4), 1608-1620. doi:10.3758/s13428-015-0672-2
 
 ### Videos
 Click below to follow through to youtube videos...

ddToolbox/CODA/CODA.m

@@ -269,9 +269,34 @@ function plotUnivariateSummaries(obj, variables, plotOptions, idNames)
 				catch
 					samples.(fieldsToGet{n}) = NaN;
 				end
+				samples.(fieldsToGet{n}) = squeeze(samples.(fieldsToGet{n}));
 			end
 		end
 
+
+		function [samples] = getSamplesAtIndex_asStochastic(obj, index, fieldsToGet)
+			samples_struct = getSamplesAtIndex_asStruct(obj, index, fieldsToGet);
+
+			% converts a structure of matricies, to a structure of
+			% Stochastics, and they are arrays of Stochastic's when
+			% appropriate
+
+			% convert samples to an array of Stochastic objects ===========
+			fieldnames = fields(samples_struct);
+			for n = 1:numel(fieldnames)
+				ncols = size(samples_struct.(fieldnames{n}),2);
+				% add samples to it
+				for col = 1:ncols
+					% create new element in object array
+					samples.(fieldnames{n})(col) = Stochastic(fieldnames{n});
+					% add samples to it
+					samples_to_add = samples_struct.(fieldnames{n})(:,col);
+					samples.(fieldnames{n})(col).addSamples( samples_to_add );
+				end
+			end
+		end
+
+
 		function [samplesMatrix] = getSamplesAtIndex_asMatrix(obj, index, fieldsToGet)
 			samples = getSamplesAtIndex_asStruct(obj, index, fieldsToGet);
 

ddToolbox/Data/Data.m

@@ -10,8 +10,14 @@
         unobservedPartipantPresent
         nExperimentFiles		% includes optional unobserved participant
 		nRealExperimentFiles	% only includes number of real experiment files
+
+		metaTable
+		% metaTable is (currently) provided by the user as an
+		% optional input into the constructor. This is meant to contain
+		% information about the overall experiment, such as participant
+		% names, id's, ages, genders, conditions etc.
 	end
-    
+	
     properties (Dependent)
         totalTrials
         groupTable % table of A, DA, B, DB, R, ID, PA, PB
@@ -30,6 +36,8 @@
 			p.addRequired('dataFolder',@isstr);
 			p.FunctionName = mfilename;
 			p.addParameter('files',[],@(x) iscellstr(x)|ischar(x));
+			p.addParameter('metaTable', table(), @istable);
+			p.addParameter('metaTableFile', [], @isstr);
 			p.parse(dataFolder, varargin{:});
 
 			assert(~isempty(p.Results.files), 'no filenames provided under ''files'' input argument')
@@ -44,6 +52,46 @@
 			if ~isempty(p.Results.files)
 				obj = obj.importAllFiles(p.Results.files);
 			end
+
+
+			%% Parse things related to metaTable ==========================
+			% Note that the first column of the file MUST be the filenames
+
+			% default behaviour if no arguments:
+			if isempty(p.Results.metaTable) && isempty(p.Results.metaTableFile)
+				obj.metaTable = table(p.Results.files', 'RowNames',p.Results.files');
+				obj.metaTable.Properties.VariableNames{1} = 'filename';
+			end
+
+			% error if both file and table supplied
+			if ~isempty(p.Results.metaTable) && ~isempty(p.Results.metaTableFile)
+				error('Pass in EITHER a table OR a filename to an table, not both.')
+			end
+
+			% user supplied metaTable
+			if ~isempty(p.Results.metaTable)
+				% Check 1: metaTable must have Row names equal to those in files
+				assert(numel(p.Results.metaTable.Properties.RowNames) == numel(p.Results.files),...
+					'metaTable must have same number of rows as number of files passed in')
+				% Check 2: they actually match up
+				warning('Implement this validation check for extra safety.')
+				obj.metaTable = p.Results.metaTable;
+			end
+
+			% user supplied path to a .csv file with experimental info
+
+			if ~isempty(p.Results.metaTableFile)
+				metaTable = readtable(p.Results.metaTableFile);
+				% set rownames equal to the column called 'filename'
+				metaTable.Properties.RowNames = metaTable.filename;
+% 				% Check 1: metaTable must have Row names equal to those in files
+% 				assert(numel(metaTable.Row) == numel(p.Results.files),...
+% 					'metaTable must have same number of rows as number of files passed in')
+% 				% Check 2: they actually match up
+% 				warning('Implement this validation check for extra safety.')
+				obj.metaTable = metaTable;
+			end
+
 		end
 
 
@@ -85,7 +133,8 @@ function exportGroupDataFileToDisk(obj)
 				error('Have already added unobserved participant')
 			end
 
-			obj.filenames{obj.nRealExperimentFiles+1} = unobservedParticipantString;
+			obj.filenames{obj.nRealExperimentFiles+1}	   = unobservedParticipantString;
+			obj.filenames_full{obj.nRealExperimentFiles+1} = unobservedParticipantString;
 			obj.participantIDs{obj.nRealExperimentFiles+1} = unobservedParticipantString;
 
 			obj.nExperimentFiles = obj.nExperimentFiles + 1;
@@ -176,16 +225,7 @@ function exportGroupDataFileToDisk(obj)
 		end
 
 		function nTrials = getTrialsForThisParticant(obj, p)
-			% TODO: really need to get a better solution that this special
-			% case nonsense for the unobserved group participant with no
-			% data
-			if p > numel(obj.experiment)
-				% asking for an experiment which doesnt exist. Probably
-				% happening because of the group-level estimate
-				nTrials = [];
-			else
-				nTrials = obj.experiment(p).getTrialsForThisParticant;
-			end
+            nTrials = height(obj.experiment(p).getDataAsTable());
 		end
 
 		function pTable = getRawDataTableForParticipant(obj, p)
@@ -206,14 +246,14 @@ function exportGroupDataFileToDisk(obj)
 			if ischar(whatIwant)
 				switch whatIwant
 					case{'all'}
-						names = obj.filenames;
+						names = obj.filenames_full;
 					case{'experiments'}
-						names = obj.filenames([1:obj.nRealExperimentFiles]);
+						names = obj.filenames_full([1:obj.nRealExperimentFiles]);
 					case{'group'}
 						if ~obj.unobservedPartipantPresent
 							error('Asking for group-level (unobserved) participant, but they do not exist')
 						end
-						names = obj.filenames(end);
+						names = obj.filenames_full(end);
 				end
 			elseif isnumeric(whatIwant)
 				% assume we want to index into IDnames
@@ -286,6 +326,15 @@ function exportGroupDataFileToDisk(obj)
 				groupTable = vertcat(groupTable, obj.experiment(n).getDataAsTable);
 			end
 		end
+
+		function metaTable = getMetaTable(obj)
+			metaTable = obj.metaTable;
+		end
+
+		function proportionDelayedOptionsChosen = getProportionDelayedOptionsChosen(obj,n)
+			responses = obj.getParticipantResponses(n);
+			proportionDelayedOptionsChosen = sum(responses)/numel(responses);
+		end
 
 	end
 

ddToolbox/Data/DataFile.m

@@ -19,16 +19,6 @@
 
 			% TODO: throw error if we have no data
 		end
-
-		function aTable = getDataAsTable(obj)
-			aTable = obj.datatable;
-			assert(istable(aTable))
-		end
-
-		function nTrials = getTrialsForThisParticant(obj)
-			% TODO: rename this function
-			nTrials = height(obj.datatable);
-		end
 
 		function obj = plot(obj, dataPlotType, timeUnits)
 			% This should be able to deal with:
@@ -40,9 +30,7 @@
 			else
 				timeUnitFunction = str2func(timeUnits);
 			end
-
-
-
+
 			% exit if we have got no data
 			if isempty(obj.datatable)
 				warning('Trying to plot, but have no data. This is probably due to this being the (group/unobserved) participant, who has no data. This is only an error if you are not getting data corresponding to a specific data file.')
@@ -118,8 +106,14 @@
 
 		end
 
+        % PUBLIC GETTERS =======================================================
 
-		function r = getDelayRange(obj)
+        function aTable = getDataAsTable(obj)
+            aTable = obj.datatable;
+            assert(istable(aTable))
+        end
+
+		function r = getUniqueDelays(obj)
 			try
 				r = unique(sort([obj.datatable.DA(:) ;obj.datatable.DB(:)]));
 			catch
@@ -128,7 +122,6 @@
 			end
 		end
 
-
 	end
 
 
@@ -141,48 +134,59 @@
 
 		function [x,y,z,markerCol,markerSize] = convertDataIntoMarkers_Homogenous(obj)
 			% FOR 2D DISCOUNT FUNCTIONS
-
-			% find unique experimental designs
-			D=[abs(obj.datatable.A), abs(obj.datatable.B), obj.datatable.DA, obj.datatable.DB];
-			[C, ia, ic] = unique(D,'rows');
-			%loop over unique designs (ic)
+			[ia, ic] = obj.calcUniqueDesigns();
+            markerSize = obj.calcMarkerSize(ic);
+            markerCol = obj.calcMarkerCol(ic);
+
 			for n=1:max(ic)
-				% binary set of which trials this design was used on
-				myset=ic==n;
-				% Size = number of times this design has been run
-				markerSize(n) = sum(myset);
-				% Colour = proportion of times that participant chose immediate
-				% for that design
-				markerCol(n) = sum(obj.datatable.R(myset)==0) ./ markerSize(n);
-
-				%x(n) = abs(p.Results.datatable.B( ia(n) )); % �B
-				x(n) = obj.datatable.DB( ia(n) ); % delay to get �B
+				x(n) = obj.datatable.DB( ia(n) ); % delay to get B
 				y(n) = abs(obj.datatable.A( ia(n) )) ./ abs(obj.datatable.B( ia(n) ));
 			end
-			z=[];
+			z = [];
 		end
 
 		function [x,y,z,markerCol,markerSize] = convertDataIntoMarkers_Heterogenous(obj)
 			% FOR 3D DISCOUNT FUNCTIONS
-
-
-			% find unique experimental designs
-			D=[abs(obj.datatable.A), abs(obj.datatable.B), obj.datatable.DA, obj.datatable.DB];
-			[C, ia, ic] = unique(D,'rows');
-			% loop over unique designs (ic)
+            [ia, ic] = obj.calcUniqueDesigns();
+            markerSize = obj.calcMarkerSize(ic);
+            markerCol = obj.calcMarkerCol(ic);
 			for n=1:max(ic)
-				% binary set of which trials this design was used on
-				myset=ic==n;
-				% markerSize = number of times this design has been run
-				markerSize(n) = sum(myset);
-				% Colour = proportion of times participant chose immediate for that design
-				markerCol(n) = sum(obj.datatable.R(myset)==0) ./ markerSize(n);
-
-				x(n) = abs(obj.datatable.B( ia(n) )); % �B
-				y(n) = obj.datatable.DB( ia(n) ); % delay to get �B
+				x(n) = abs(obj.datatable.B( ia(n) )); % B
+				y(n) = obj.datatable.DB( ia(n) ); % delay to get B
 				z(n) = abs(obj.datatable.A( ia(n) )) ./ abs(obj.datatable.B( ia(n) ));
 			end
 		end
-
+
+        function [ia, ic] = calcUniqueDesigns(obj)
+            designs = [abs(obj.datatable.A), abs(obj.datatable.B), obj.datatable.DA, obj.datatable.DB];
+            [C, ia, ic] = unique(designs,'rows');
+        end
+
+        function markerSize = calcMarkerSize(obj, ic)
+            for n=1:max(ic)
+                markerSize(n) = obj.nDesignOccurrences(ic, n);
+            end
+        end
+
+        function markerCol = calcMarkerCol(obj, ic)
+            % Colour = prop. of times they chose immediate, for that design
+            for n=1:max(ic)
+                markerCol(n) = obj.nChoseImmediate(obj.datatable.R, ic, n)...
+					./ obj.nDesignOccurrences(ic, n);
+            end
+		end
+
 	end
+
+    methods(Static)
+
+		function nImmediateChoicesMade = nChoseImmediate(responses, ic, n)
+			nImmediateChoicesMade = sum(responses(ic==n)==0);
+		end
+
+        function nOccurrences = nDesignOccurrences(ic, n)
+            nOccurrences = sum(ic==n);
+        end
+
+    end
 end