# Subject-Specific Scalar-on-Image Regression (S3IR)

## Overview
S3IR is a statistical method for analyzing heterogeneous neuroimaging data by incorporating subject-specific spatial patterns while maintaining population-level inference capabilities.

## Installation
1. Extract the package to your desired location
2. Run `init` in MATLAB to set up paths

## Quick Start
```matlab
% Initialize package
init;

% Run simulation example
simulation_example;

% For real data analysis
help real_data_analysis
```

## Package Structure
```
s3ir/
├── init.m                          % Package initialization
├── README.md                       % This file
├── core/                          % Core algorithm components
│   ├── S3IR.m                     % Main algorithm
│   ├── B_update.m                 % Population mask update
│   ├── tau_update.m               % Potts parameter update
│   └── region_similarity.m        % Similarity computation
├── rkhs_utils/                    % RKHS regression utilities
│   ├── gaussian_kernel_sparse.m   % Kernel matrix construction
│   ├── rkhs_regression.m          % Kernel regression
│   └── nbr_radio_based.m          % Neighborhood coordinates
├── spatial_utils/                 % Spatial processing utilities
│   ├── head_file_for_mask.m       % Neighborhood structure
│   └── potts_likelihood.m         % Potts derivatives
├── prediction_utils/              % Prediction utilities
│   ├── predict_standard.m         % Standard prediction
│   ├── predict_sensitivity.m      % Sensitivity analysis
│   └── weight_assign.m            % Pattern matching
├── sim_utils/                     % Simulation utilities
│   ├── generate_synthetic_data.m  % Data generation
│   └── random_region_assignment.m % Region assignments
└── examples/                      % Example scripts
    ├── simulation_example.m       % Simulation demo
    └── real_data_analysis.m       % Real data analysis
```

## Main Components

### Core Algorithm: `S3IR`
The main estimation algorithm implementing the EM-type procedure for fitting subject-specific scalar-on-image regression models.

### Prediction Utilities
- `predict_standard`: Standard prediction using pattern matching
- `predict_sensitivity`: Prediction with sensitivity parameter for robustness analysis

### Simulation Utilities
- `generate_synthetic_data`: Creates synthetic neuroimaging data for simulation studies
- `random_region_assignment`: Generates random region assignments

### RKHS Utilities
- `gaussian_kernel_sparse`: Constructs sparse Gaussian kernel matrices
- `rkhs_regression`: Performs kernel ridge regression with BIC selection
- `nbr_radio_based`: Generates neighborhood coordinate offsets

### Spatial Utilities
- `head_file_for_mask`: Creates spatial neighborhood structures
- `potts_likelihood`: Computes Potts model derivatives

## Key Features
- Handles subject heterogeneity in neuroimaging data
- Spatial regularization using Potts model
- RKHS-based coefficient estimation with Gaussian kernels
- Automatic hyperparameter selection via BIC
- Pattern matching prediction framework
- Support for 1D, 2D, and 3D images

## Usage Examples

### Basic Simulation
```matlab
% Set up simulation parameters
options.itermax1 = 50;
options.itermax2 = 50;
options.lambda_sparse = 100;
options.tau = 10;
options.RKHS_band_width_s = 1:2;
options.RKHS_lambda_s = 2.^(-5:1:5);

% Generate synthetic data
[data, ~, ~, ~, ~] = generate_synthetic_data(200, 0.5, [48, 48], 3);

% Fit S3IR model
model = S3IR(data, options);

% Make predictions
[y_pred, assignments] = predict_standard(model, data.X_new);
```

### Real Data Analysis
```matlab
% Load your neuroimaging data
% data.X: [n x q] covariate matrix
% data.Y: [n x 1] response vector  
% data.imagesize: image dimensions
% data.mask: valid pixel mask

% Set algorithm options
options = struct();
options.itermax1 = 30;
options.lambda_sparse = 10;
options.tau = 10;
options.RKHS_band_width_s = 1:2;
options.RKHS_lambda_s = 2.^(-5:1:5);

% Fit model
model = S3IR(data, options);

% Evaluate performance
predictions = predict_standard(model, test_data.X);
RMSPE = sqrt(mean((predictions - test_data.Y).^2));
correlation = corr(predictions, test_data.Y);
```

## Real Data Requirements

### ADNI Data Access
The real data analysis example (`real_data_analysis.m`) is designed for ADNI (Alzheimer's Disease Neuroimaging Initiative) data, but the actual data cannot be provided due to privacy restrictions. To use this functionality:

1. **Obtain ADNI Access**: Apply for data access at http://adni.loni.usc.edu/
   - Registration and approval required
   - Institutional agreements may be necessary
   - Processing time typically 1-2 weeks

2. **Required Data Files**:
   - `ANTS_shrinked_registered.mat`: Preprocessed MRI data structure containing:
     - `ANTS.imdata`: [749 × q] imaging data matrix (RAVENS maps)
     - `ANTS.mask`: [128×128×128] brain mask indicating valid voxels
     - `ANTS.table`: Clinical assessments with variables:
       - `DX_bl`: Baseline diagnosis (Normal Control, MCI, Dementia)
       - `MMSE_bl`, `MMSE_06`: MMSE scores at baseline and 6-month follow-up
       - `ADAS11_bl`, `ADAS11_06`: ADAS-Cog 11-item scores
       - `ADAS13_bl`, `ADAS13_06`: ADAS-Cog 13-item scores
   - `ANTS_tr.txt`: Training set indices for cross-validation folds
   - `ANTS_te.txt`: Test set indices for cross-validation folds

3. **Data Preprocessing Pipeline**:
   The ADNI data should be preprocessed using ANTs (Advanced Normalization Tools):
   - Anterior commissure and posterior commissure correction
   - N2 bias field correction  
   - Skull-stripping
   - Intensity inhomogeneity correction
   - Cerebellum removal
   - Segmentation and registration
   - RAVENS map generation (tissue density quantification)

4. **Study Population**:
   The analysis expects 749 participants including:
   - 165 Alzheimer's Disease patients
   - 366 Mild Cognitive Impairment patients  
   - 218 Healthy controls

### Alternative Data Sources
If ADNI data is not available, the framework can be adapted for other neuroimaging datasets by:
- Ensuring similar preprocessing pipeline
- Adapting the data structure format in `real_data_analysis.m`
- Modifying response variable selections as needed

### Data Format Requirements
Your data structure should match:
```matlab
data.imagesize = [D1, D2, D3];  % Image dimensions
data.mask = logical_mask;       % Valid voxel mask
data.X = subject_images;        % [n_subjects × n_voxels] imaging data
data.Y = clinical_scores;       % [n_subjects × 1] outcome variable
```

## Requirements
- MATLAB R2012b or later
- Image Processing Toolbox (recommended)
- Statistics and Machine Learning Toolbox (recommended)

## Algorithm Overview
S3IR addresses heterogeneity in neuroimaging data through:

1. **Population-level modeling**: Estimates a common disease map using RKHS regression
2. **Spatial regularization**: Uses Potts model to enforce spatial smoothness
3. **Subject-specific adaptation**: Binary masking images capture individual heterogeneity
4. **Pattern matching prediction**: Determines active regions for new subjects

## Model Parameters

### Core Parameters
- `itermax1`: Maximum outer EM iterations (default: 50)
- `itermax2`: Maximum inner iterations for mask update (default: 50)
- `lambda_sparse`: Sparsity penalty for region assignment (default: 100)
- `tau`: Potts model spatial smoothness parameter (default: 10)

### RKHS Parameters
- `RKHS_band_width_s`: Gaussian kernel bandwidth candidates
- `RKHS_lambda_s`: Ridge penalty candidates for kernel regression

### Region Detection
- `num_region`: Maximum number of regions to detect (default: 10)
- `region_size_thres`: Minimum region size threshold (default: 13)

## Performance Metrics
The package provides several evaluation metrics:
- **REE**: Relative Estimation Error for coefficient recovery
- **RMSPE**: Root Mean Squared Prediction Error
- **Correlation**: Pearson correlation between predicted and observed responses
- **AIC/BIC**: Model selection criteria

## Troubleshooting

### Common Issues
1. **Memory errors**: Reduce image size or use sparser kernel matrices
2. **Convergence problems**: Adjust iteration limits or regularization parameters
3. **No regions detected**: Check mask validity and reduce region size threshold

### Performance Tips
- Use appropriate image preprocessing (skull stripping, normalization)
- Start with conservative parameter settings for real data
- Monitor convergence through iteration output messages

## Citation
Please cite the original research when using this code.

## License
This software is provided for research purposes. Please refer to the original publication for usage terms.