博士学位論文

Paulus Kapundja Shigwedha

氏 名 (本 籍 )
学位の種類

博

士（工学）

学位記番号

甲第１７０号

学位授与年月日

令和

3

年

4

学位授与の要件

学位規則第４条第

月

26

日

1

項

学位論文題目
Methods to reduce the bias and enhance
the contrast of receptor density inposi tron emission tomography imaging

論文審査委員

主査

福岡

豊

副査

於保

英作

〃

斎藤

秀俊

〃

木村

裕一

〃
〃
〃

工学院大学大学院

© Paulus Kapundja SHIGWEDHA

2020

i

Dedications:
I dedicate this thesis to my late father, Tomas Shikongeni SHIGWEDHA, and to my two younger
brothers, Leevi Panduleni SHIGWEDHA and Esra Gandja-Ayihe SHIGWEDHA.

ii

Acknowledgements
Firstly, I would like to show my gratitude towards the Japanese government, which, through its
Ministry of Education, Culture, Sports, Science and Technology (MEXT), offered me a scholarship to come further my studies in this wonderful nation of the rising sun, Japan.

I am also immensely grateful to my supervisor [Yutaka FUKUOKA, Professor, Faculty of Engineering, Department of Electrical and Electronic Engineering, Kogakuin University of Technology and Engineering (KUTE), Tokyo, Japan] for his vast guidance and motivation during
my studies. He has been continuously supportive, and that made it possible for me to achieve
the goal of my research. I appreciate the input by our external supervisor [Yuichi KIMURA,
Professor, Department of Computational Systems Biology, Faculty of Biology-Oriented Scienceand Technology, Kindai University, Wakayama, Japan], he provided me with his well established
MATLAB® library for the LGA execution, and his valuable assistance throughout the course ofthis research is fully acknowledged. In addition, I appreciate the members of the Thesis Committee, co-supervisor [Eisaku OHO, Professor, Faculty of Engineering, Department of Electricaland Electronic Engineering, KUTE, Tokyo, Japan] and [Hidetoshi SAITO, Professor, Facultyof Engineering, Department of Electrical and Electronic Engineering, KUTE, Tokyo, Japan],for their valuable time as well comments that greatly improved the content of this thesis. I
am also thankful to my fellow labmates, in the Laboratory for Biomedical Informatics (LBMI),
KUTE, academic years, 2015–2020, for their help and discussions in the laboratory. Further, I
would also like to appreciate the support provided by KUTE’s Student Support Division, andthe entire University community throughout my studies.

My family has always been there for me during my studies. I would like to thank my motherand my late father, and my entire family and friends back home for their unlimited love, guidance and support in all forms. I would also like to thank the Namibian personnel in Japan for
their productive moments we have had during my studies. I am as well thankful to my friends
in Japan, especially Mr. Hiroki HARA for the time he has spent helping me in various ways,
especially with social situations during the course of my stay in Japan.

I reserved special thanks for my best friend Hilda Theofilus. She has been encouraging andalways ready to listen and attend to me whenever I needed her.
iii

Contents
1 Introduction

1

2 Positron emission tomography, compartmental model and the Logan graphical analysis

5

2.1

Positron emission tomography . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

2.1.1

Positron emission tomography detection system . . . . . . . . . . . . . . .

7

2.1.2

Compartmental model and the Logan graphical analysis . . . . . . . . . . 10

2.2

Logan graphical analysis using a reference region . . . . . . . . . . . . . . . . . . 19

3 Reducing the bias in the BPN D estimates using least-squares cubic linear
regression
3.1

3.2

31

Least-squares cubic linear regression . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.1.1

Simulation studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.1.2

Clinical 11 C-PiB PET data studies . . . . . . . . . . . . . . . . . . . . . . 50

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4 Tissue time-activity curves denoising techniques for improving image contrast
4.1

4.2

4.3

58

Principal component analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.1.1

Simulation studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

4.1.2

Clinical 11 C-PiB PET data studies . . . . . . . . . . . . . . . . . . . . . . 68

Correlated component analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.2.1

Simulation studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4.2.2

Clinical 11 C-PiB PET data studies . . . . . . . . . . . . . . . . . . . . . . 83

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

iv

PK SHIGWEDHACONTENTS

5 Summary and future directions

93

A Reducing the least-squares cubic equation to commonly known regression
methods under specific assumptions

95

A.1 Only Yi is subject to errors (i.e., Xi is error-free): . . . . . . . . . . . . . . . . . . 96
A.2 Only Xi is subject to errors (i.e., Yi is error-free): . . . . . . . . . . . . . . . . . . 96
A.3 Errors present in both X and Y variables: . . . . . . . . . . . . . . . . . . . . . . 96
B Noise model98

C Contrast evaluation of simulation data

99

Bibliography

102

工学院大学

v

KOGAKUIN UNIVERSITY

Chapter 1

Introduction
Positron emission tomography (PET) is a physiological imaging technique in nuclear medicine.
PET is based on the fundamentals of nuclear dacay, specifically beta-plus (β + ) decay. Beta-plus
decay is a nuclear disintegration process that occurs in a nucleus with more number of protons
than neutrons. In a β + decay, a proton is transformed into a neutron, and then a positron anda neutrino are simultaneously released [1]. PET is used to evaluate physiological processes of
biological systems, which is applicable in diagnosing, staging and treatment of various diseases,
such as various heart and brain diseases and cancers. In a PET study, a β + radioactive substance, referred to as a radiotracer, is administered in vivo. When the administered radiotracer
undergoes β + decay in vivo, emitted positrons annihilates with the electrons in tissues. The
positron-electron annihilation gives off two γ-rays, which are approximately in opposite directions. The γ-rays are then detected to construct the concentration of the radiotracer in tissues
over time; these data are referred to as tissue time-activity curves (tTACs). Various methods
such as the Logan graphical analysis (LGA) can then be applied to tTACs for quantitative
analysis of the radiotracer kinetics.

LGA can be thought of as a linearization technique for PET data. LGA transforms PET tTACs
into linear relationships [2, 3]. Physiologically interpretable quantities, such as distribution
volume (DV ), distribution volume ratio (DV R) and the non-displaceable binding potential(BPN D ) can be obtained from the slope of the LGA linear relationship. LGA is computationallyefficient, in that it is easy to implement and fast to compute, owing to which it has been wellaccepted amongst researchers.

1

PK SHIGWEDHAINTRODUCTION

Problem statement and motivation
Traditionally, the LGA slope is estimated by the common ordinary least-squares (OLS) regression. However, in the presence of noise in PET data, it has been shown that the OLS-based
estimates of the LGA slope are underestimated [4, 5, 6, 7]. In addition, this negative bias has
been observed to increase with the magnitude of both noise and the LGA slope [4, 5, 6, 7].
The underestimation is due to that both the LGA variables are contaminated with a correlated
noise, but OLS only accounts for the noise in the response variable. Specifically, the noise in the
LGA variables is attributed to the noisy term, C(t) (see Eqs. (2.21, 2.40)), in the denominatorsof both of the two LGA variables. Accordingly, the noise due to C(t) propagates correlatively
in the LGA variables, presenting a correlated errors-in-variables (EIV) problem. This explicitly
makes OLS unsuitable for estimating the LGA slope.

Various methods have been proposed to reduce the influence of the noise in PET data on theestimated slope. These methods range from tTAC smoothing methods, such as the generalized
linear least-squares (GLLS) [8] technique, to methods based on a rearrangement of the original
LGA equation, such as the multilinear reference tissue models (MRTM 1 & 2) [9, 10]. Other
methods including orthogonal distance regression (ODR) [11] simply employ a different line
estimation technique. These methods, however, have been found to either only slightly reducethe bias or reduce the bias whilst causing variations in the estimates [5, 6]. This necessitates
the need for further studies. This study aims to establish a method to reduce the bias in theestimates of the LGA slope, and improve the contrast of the resulting parametric images. In[11], it has been demonstrated that the regression method used to estimate the Logan slope
influences the resulting bias. On this basis, this study employs an alternative linear regression
method referred to as least-squares cubic (LSC) to estimate the LGA slope. LSC accounts fornoise in both variables, by minimizing the squared residuals in both the predictor and response
variables [12, 13, 14]. In addition, LSC also incorporates the correlation of the noise in the variables. To further improve the LSC estimates, LSC is for the first time combined with tTACs
denoising techniques, principal component analysis (PCA) and correlated component analysis(CorrCA).

PCA is a dimensionality reduction technique used mostly for feature extraction and noise filtering. PCA transforms variables into new sets of variables that are linear functions of the original
工学院大学

2

KOGAKUIN UNIVERSITYPK SHIGWEDHAINTRODUCTION

variables [15]. The new sets of variables are uncorrelated and are defined by sets of orthogonalbasis vectors and principal components (PCs) that optimally describe the variance in the data.
The PCs and the corresponding variance are obtained by solving for the eigenvectors and eigenvalues of the covariance matrix of the data. The eigenvectors are the orthogonal basis vectors,
which define the PCs, and the eigenvalues are the corresponding variances.

PCA is a common tool in PET parametric imaging; it has been previously used to denoisetTACs in [16], and it showed minimal effects on the variance. Therefore, this effect is employedin this study. PCA has also been used for clustering in [17]. Further, PCA has also been used
for noise filtering in PET studies with irreversible radiotracers [18]. Here, LSC and PCA areused together for the first time, and the effect of this fusion on the image quality was assessed
not only in the usual manner of visual inspection of the images but also numerically in termsof the contrast between grey matter and white matter. In other words, this study utilizes the
combined effects of LSC’s minimal bias and PCA’s precision.

Like PCA, CorrCA is also a feature extraction method. CorrCA operates by identifying components that are maximally correlated between repetitions in multivariate data [19]. Specifically,
CorrCA maximizes the ratio of between-repetition to within-repetition covariance [19]. In thecontext of our application to PET data, our repetitions are along the slices dimension. CorrCA
will therefore maximize the ratio of between-slices to within-slices covariance. This ratio is referred to as inter-subject correlation (ISC). Specifically in our context, it translates to inter-slice
correlation.

Maximizing the between-slices to within-slices covariances maximizes the mean-over-varianceacross slices, which has been asserted to define a signal-to-noise ratio measure [19]. CorrCA
capitalize on the ability to simultaneously operate ”within individual slices” and ”across all
slices”, without arranging the slices in one plane. That means operating directly on the 3Darray, unlike PCA which requires the slices to be arranged side-by-side into a 2D-array. This
could be advantageous over PCA, which means CorrCA can possibly provide a better noise
filter method in comparison to PCA. This study therefore also seek to introduce CorrCA to
PET parametric imaging and assess the performance of CorrCA in comparison to PCA. Hereafter, the fusions of LSC with PCA and CorrCA are respectively denoted as LSC–PCA and

工学院大学

3

KOGAKUIN UNIVERSITYPK SHIGWEDHAINTRODUCTION

LSC–CorrCA. Similarly, using OLS with PCA and CorrCA are herein referred to as OLS–PCA
and OLS–CorrCA, respectively.

This study was conducted on both simulation and real data. Two sets of simulation data wereused, one mimicking 11 C-carfentanil (CFN), and the other mimicking 11 C Pittsburgh compound
B (11 C-PiB). The real data are for 11 C-PiB PET brain images. Both CFN and 11 C-PiB havespecific binding sites in the brain. CFN binds to µ-opioid receptors, and 11 C-PiB binds to beta
amyloid (Aβ) plaques.

Thesis outline
The next chapter (Chapter 2) goes into more detail about PET. In the subsections, the compartmental model, derivation of the LGA from the compartment model, as well as the interpretationof the parameters estimated from the LGA are discussed. OLS, ODR and MRTM2 are also
discussed. Chapter 3 gives the derivation and analysis of LSC. The analysis of LSC is also
given specifically in terms of the LGA variables. The results of the study of LSC-based LGA incomparison to OLS-based LGA and MRTM2 are presented in this chapter (Chapter 3). These
are the results reported in the published paper: [20]. In Chapter 4, the mathematical detailsof both PCA and CorrCA are presented. The results of the study of LSC–CorrCA, LSC–PCA,LSC and OLS are given in this chapter (Chapter 4). Chapter 5 summarizes this thesis anddiscusses future directions.

工学院大学

4

KOGAKUIN UNIVERSITY

Chapter 2

Positron emission tomography,
compartmental model and the Logangraphical analysis
2.1

Positron emission tomography

PET is referred to as a functional imaging method. This is to differentiate it from other imaging
techniques such as magnetic resonance imaging (MRI) which are defined as structural imaging.
A functional imaging helps image the functioning status of cells in vivo, whereas structural
imaging is used to examine the structural shapes of tissues and organs. PET is used to visualize
physiological functioning of tissue cells by monitoring the interaction of the radiotracer withthe system’s biochemical processes and activities such as glucose metabolism, oxygen consumption and blood flow. PET is advantageous over other imaging techniques in that it can detect
diseases from their early onsets. In a similar way PET can also detect immediate response of apatient to treatments.

A radiotracer is made to have an affinity to a particular region of interest (ROI). Therefore,
eventually after administration, the radiotracer will accumulate in the ROI. As the radionuclide
undergoes β + decay in the ROI, the time-course of radioactivity in tissues is recorded to be used
to study the physiological activities of interest in the ROI. Figure 2.1 shows a simple schemeof PET imaging. Radiotracers are made of two components; a radionuclide and a biologicalmolecule. This forms a biological radioactive compound. The biological molecule describes the5

PK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGAaffinity aspects of the radiotracer.

Equation (2.1) [1] represents a β + decay process. The parent nucleus A
Z XN has more numberof protons than neutrons, i.e., Z > (A − Z = N ). A proton in the proton-rich parent nucleus
is transformed into a neutron, accompanied by a simultaneous emission of a positron (e+ ) anda neutrino (v). The resulting daughter nucleus, A
Z−1 YN +1 , has ”one-less number of protons”
and ”one-more number of neutrons” than the parent nucleus. The mass number (A) of thedaughter nuclide remains the same as that of the parent nuclide, i.e., the daughter nuclide is an
isobar of the parent nuclide —different chemical elements with same number of nucleons (A)and different number of protons. The daughter nuclide can be a stable isotope or a radioactive
isotope which is closer (than the parent nuclide) to the stability region on the chart of nuclides.
AA+
Z XN → Z−1 YN +1 + e + v.

(2.1)

The β + decay equation for 11 C is represented in Eq. (2.2). 11 C is the radionuclide label for
11

C-PiB radiotracer, which is utilized for clinical 11 C-PiB PET images analyzed in this study.
11
11+
6 C5 → 5 B6 + e + v + 0.96 MeV.

(2.2)

Equation (2.2) means that the unstable 11 C radionuclide undergoes β + decay to form the stable
isotope 11 B, and a total of 0.96 MeV is released —which is mostly carried away as kinetic energyof the released particles. 11 C has a half-life of 20.364 minutes.

Beta-plus decay is in contrast to another form of beta decay, referred to as beta-minus (β − )decay, in which a neutron is transformed into a proton, accompanied by a simultaneous emissionof an electron (e− ) and an anti-neutrino (v̄). Beta-minus decay can be represented by Eq. (2.3)[1] below:AA−
Z XN → Z+1 YN −1 + e + v̄.

(2.3)

With reference to the type of beta particles released in β + and β − decay, they can be referred toas positron emission and electron emission, respectively. The two beta decay processes canbe further explained interms of “up quarks” and “down quarks” —which are the fundamental
particles that forms the nucleons.

PET radiotracers are commonly made of low elements short–lived isotopes such as 11 C, 18 F,
15O and 13 N. Due to its effectiveness in cancer detection and staging, [18 F]-fluorodeoxyglucose

工学院大学

6

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA(FDG) is the widely used radiotracer [1, 21]. Depending on factors such as the ROI, the radiotracer being used and the demographics of the subject, a radiotracer can be administeredby either injection, swallowing or inhalation. The main purpose of a radiotracer is to trace thephysiological process of interest in the ROI. Therefore, it is administered in such an amountthat it does not perturb the physiological processes it is being used trace [22, 21].

Being a glucose analog, [18 F]-FDG, is transported into cells by the same mechanism as glucose.
Intracellular, glucose is phosphorydated into glycolysis. Unlike glucose, [18 F]-FDG cannot enter glycolysis. This is because [18 F]-FDG does not have the 2-hydroxyl group (-OH) requiredfor glycolysis [23, 24]. [18 F]-FDG therefore gets trapped in cells in the form of [18 F]-FDG-6Phosphate [23]. [18 F]-FDG-6-Phosphate exits the cells via β + radioactive decay. Imaging of the
distribution of the [18 F]-FDG by PET can therefore be used to trace the distribution of glucose
uptake and phosphorylation.

The clinical data used in this study are of 11 C-PiB radiotracer. 11 C-PiB binds to Aβ plaques inthe brain [25, 26], which are associated with Alzheimer’s disease (AD) and and mild cognitive
impairment (MCI). Quantification of 11 C-PiB binding therefore allows for visualization andquantification of in vivo Aβ deposition.

2.1.1

Positron emission tomography detection system

PET data acquisition is based on measuring the total energy lost by annihilation photons as
they interact with the detector equipment. Annihilation photon detectors form a ring surrounding the subject to undergo PET imaging. The detection ring depicted in Fig. 2.1 is made of
scintillation materials mounted onto photomultipliers tubes (PMTs) [27]. Photon detection isthen based on their interaction with the scintillator material. The scintillator material is a
specially-designed crystal which emits visible light upon interaction with annihilation photons[1]. This process is referred to as scintillation. In short, in the best case scenario, when there isno photoelectric or scatter events, annihilation photons reach and directly excite the scintillator
valence electrons. The excited electrons get into the conduction band, and return (de-excite)
to the ground state by releasing some energy via photon emission [1]. Scintillation materials
are modified by adding impurities to pure crystals such that the scintillation photons are in thevisible electromagnetic spectrum, i.e., luminescence [1]. The PMTs then detect the scintillation

工学院大学

7

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

'

$

&

%

Figure 2.1: A simple PET schematic diagram.

photons by ionization mechanism of photoelectric effect. Two photons detected within a presetcoincidence time window (usually a few nanoseconds (1 – 12) ) are recorded and used to determine a line of response (the line connecting the centers of the two detectors.) [1, 27, 28]. The
number of scintillation photons produced in the scintillator medium, and the subsequent electriccharge produced in the PMT is proportional to the energy deposited by the annihilation photon
in the scintillator medium. These data can then be used to evaluate radioactivity concentrationin tissues.

The four main properties of scintillator materials are, stopping power for annihilation photon(511 keV), signal decay time, light output, and the intrinsic energy resolution [1]. The stopping
power is defined by the inverse of the mean distance a photon can travel before depositing its
energy in the detector material [1, 28]. This distance is termed attenuation length [1], and it isdependent on the density and effective atomic number of the scintillator material. The high thenumber of photons depositing their energy in the scintillator material the better, i.e., a scintillator material with a short attenuation length is preferable. The decay time of a scintillator

工学院大学

8

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

describes the rate at which the valence electrons excite and de-excite following annihilation photon absorption, i.e., the rate of photoluminescence. During this period the detector is referred
to as being dead since it cannot process all other coincidence events [28], and are therefore lost
—this time period is referred to as dead time. Decay time is usually in the range of hundredsof nanoseconds (e.g., 230 nanoseconds for Sodium Iodide doped with Thallium (NaI(Tl)) scintillator crystal [28]). The manner in which the decay time is associated with the efficiency ofthe scintillator material is that, short decay time will allow discrete cascades of electrons tobe processed at high count rates, and minimize the number of random coincidences [1]. Light
output refers to the number of scintillation photons produced by an incident annihilation photon. High light output allows for a good spatial resolution and good energy resolution [1, 28].
High energy resolution allows for improved discrimination of scatter and true events [1, 28].
The energy resolution of a scintillator is also dependent on the intrinsic energy resolution of thescintillator material, which is dependent on the inhomogeneities of the crystal and its associatedproperties, such as light-output uniformity [1].

There are many types of scintillator materials used in PET detection systems. One such is the
NaI(Tl) crystal mentioned earlier. The NaI(Tl) scintillator has a high light output, providing
a good energy and spatial resolution [1, 28]. However, NaI(Tl) scintillator crystal has a slow
decay time, which means increased detector dead time, as well as increased random coincidences[1]. Because of its low density, NaI(Tl) also has a low stopping power than the other scintillator
materials [1, 28]. One other scintillator material is the Gadolinium Oxyorthosilicate doped with
Cerium (GSO). GSO has a lower stopping power and energy output, but it is mostly emplyedfor its improved energy resolution and uniform light output in comparison to other scintillator
materials [1]. Another scintillator material is Barium Flouride (BaF2 ). BaF2 is noted for itsvery short decay time (600 picoseconds), an excellent feature that makes up for the effects of
its low stopping power [1].

PMTs are comprised of an enclosed vacuum incorporated with a photo-cathode at the entrance
window [1]. The vacuum tube of the PMT consists of dynodes, arranged at consecutively
increasing potentials. An incident scintillation photon interacts with the atoms of the photocathode layer by photoelectric effect [1]. The resulting photo-electron is accelerated to the firstdynode by an applied electric field [1]. This accelerated photo-electron will cause further mul-

工学院大学

9

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

tiple electron emissions at the dynode. This trail of electron acceleration and further multiple
secondary emissions continues at consecutive dynodes, leading to an amplified electric signal.
The final component of the dynode series is an anode through which the collected electric charge
leaves the vacuum tube for signal processing.

In PET data acquisition, three main kinds of possible coincidence events, —namely, true, random, and scatter coincidences— can be recorded. True coincidences are the events of interestin a PET study. These events occur when two photons coming from one annihilation event aredetected by two oppositely located detectors, without encountering any major interactions withthe tissues in the field of view. Random coincidences are accidental detection of two unrelated
photons as coincidence photons, detected within the preset coincidence detection window. Twophotons from a single annihilation event that are detected after either one or both have undergone a major scatter event results in what is termed scatter coincidences. Random and scatter
events are some of the major contribution to sources of errors in PET data, in addition tofactors such as positron range. Positron range refers to the distance traveled by a positron after
emission from the nucleus before annihilation [1]. The aim of PET measurement is to determine
the position of the nucleus in which beta-plus decay took place, therefore if the annihilation
event took place at some distance from the nucleus then that will add to errors in PET data.

2.1.2

Compartmental model and the Logan graphical analysis

To better understand the LGA, one needs to understand it in terms of its mathematical construction. LGA is derived from a well known concept in kinetic modeling referred to as compartmental model. Compartmental model forms the foundation of PET kinetic analysis. Speciffically,
compartmental model is used to describe the kinetics of a radiotracer in vivo via a set of differential equations. In compartmental model, it is assumed that there are discrete physiologicalregions in which the radiotracer is concentrated [21, 29]. These physiological regions are referred
to as compartments. Applying differential Calculus, one obtains a set of differential equations
describing the rates of change in the concentration of the radiotracer in the compartments withtime. The LGA is derived from the differential equations of the compartmental model, giving
a linear relationship of which the slope is a physiologically interpretable quantity.

Figure 2.2 shows schematic figures for compartmental models; Two compartments (one tissue)

工学院大学

10

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

model in 2.2(a); Three compartments (two tissues) model in 2.2(b); and Four compartments(three tissues) model in 2.2(c). The compartments CP , CF , CN , and CS represent the concentration of the radiotracer in arterial plasma, bound-free radiotracer, non-specifically boundradiotracer, and specifically bound radiotracer, respectively. The naming of the compartments
are self-explaining. The arterial plasma compartment refers to the radiotracer in the arterial
plasma, which has not yet entered the tissues. Specifically bound compartment refers to theradiotracer specifically bound to target receptors in tissues. Non-specifically bound refers to
the radiotracer that is bound to receptors other than the target receptors.

The transportation of the radiotracer between neighboring compartments is described by theexchange rate constants, K1 [ml g −1 min−1 ] and k2–4 [min−1 ] [30, 21]. K1 denotes the delivery
rate of the radiotracer from the arterial plasma to tissues; k2 is the clearance rate from thetissues back to the plasma; k3 and k4 are the respective association and dissociation rates of theradiotracer to and from the specific binding sites; and k5 and k6 are the respective associationand dissociation rates to and from the non-specific binding sites. The pairs of the exchange rate
constants and the difference of the radiotracer concentration of the two involved compartments
are assumed to be linearly related [21]. This relationship makes it possible to form the differential equations.

The concentrations of the radiotracer in the tissue compartments at time t [min] for the onetissue compartment model in Fig. 2.2(a) are described by the differential equation below,
dCt (t)= K1 CP (t) − k2 Ct (t).
dt

(2.4)

Similarly, the concentrations of the radiotracer in the tissue compartments for a two- (Fig.
2.2(b)) and three- (Fig. 2.2(c)) tissues models are respectively described by,
dCF (t)= K1 CP (t) + k4 CS (t) − (k2 + k3 )CF (t)dt(2.5)
dCS (t)= k3 CF (t) − k4 CS (t),dt
and,dCF (t)= K1 CP (t) + k6 CN (t) + k4 CS (t) − (k2 + k3 + k5 )CF (t)dtdCN (t)= k5 CF (t) − k6 CN (t)
dt

(2.6)

dCS (t)= k3 CF (t) − k4 CS (t).
dt
工学院大学

11

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA(a) Two compartments (one-tissue) model.

(b) Three compartments (two-tissues) model.

(c) Four compartments (three-tissues) model.

Figure 2.2: Schematic figures of the (a) one-, (b) two- and (c) three-tissues compartmental
models.

工学院大学

12

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

CP (t), CF (t), CS (t) and CN (t) denote the concentrations of radioactivity in the respective compartments at time t.

PET data are not recorded separately for the compartments. The recorded PET data represents
a sum of the radioactivity in all tissue compartments. Denoting the observed PET data by
CP ET (t), for Fig. 2.2(c) we can write,
CP ET (t) = CF (t) + CS (t) + CN (t).

(2.7)

Non-linear least-squares method can be used to fit measured tTACs to the compartmental
model to estimate the kinetic parameters (exchange rate constants) [21], using arterial plasma
time-activity curve, CP (t), as input function. Estimated kinetic parameters are then used toevaluate physiological processes of the target tissues [21]. The non-linear approach is dependent
on the number of components. Graphical analysis methods such as the LGA provides a great
advantage because they are independent of the assumed model (i.e. number of compartments)[6, 21, 2]. This makes it possible to be able to obtain parameters such as DV, DVR and BPN D
irrespective of the assumed compartmental model.

The procedure below demonstrate the linearization of the compartmental model equations to
obtain the LGA. The LGA forms a simple linear relationship of two variables, and it is the slopeof this relationship, which is a function of the kinetic parameters, that defines the physiological
functioning of the target tissues.

The differential equations in (2.4), (2.5) and (2.6) can be represented in a general form as,dA(t)= KA(t) + CP (t)Q.
dt

(2.8)

A represents a column vector of the radioactivity concentrations in tissue compartments, K
represents a matrix made of the kinetic parameters, k2−6 , and Q denotes K1 . Writing these

工学院大学

13

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

expressions, A, K and, Q, out; for the one-tissue model we get,
A = Ct (t)(2.9)

K = −k2
Q = K1 .
The two-tissues model takes the form,A=h

CF (t)CS (t)−(k2 + k3 )

K =

k4

h


(2.10)


−k4

k3
Q=

iT

K1

0

iT

,

and similarly, the three-tissues model has the form,A=h

CF (t)CS (t)CN (t)

iT




−(k
+
k
+
k
)
k
k
2
3
5
4
6 



K =
k3
−k4
0 




k5
0
−k6
Q=

h

K1

0

0

iT

(2.11).

The .T represents matrix transpose.

The total radioactivity concentration in tissues, Ct (t), is given by,Ct (t),
for a one-tissue modelCt (t) = CP ET (t) = CF (t) + CS (t),for a two-tissues model ,






CF (t) + CS (t) + CN (t) for a three-tissues model(2.12)

such that the total radioactivity concentration in the whole ROI is,
CT (t) = Ct (t) + VP CP (t),(2.13)

where VP denotes the effective plasma in tissues.

工学院大学

14

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

Let us define a unit vector U as,




1
for a one-tissue model





 






 1 


 


  for a two-tissues model 1
U =  
,






 1 


 


 


 1  for a three-tissues model


 


 


 


 1
so that we can write Ct (t) as
Ct (t) = U T A(t).

(2.14)

dA(t) = K A(t) dt + CP (t)Q dt.

(2.15)

Equation (2.8) can be rearranged into,

Re-arranging Eq. (2.15) and integrating gives,
A(t)dt = −K −1 QCP (t)dt + K −1 dA(t)Z tA(u)du = −K

−1

(2.16)Z t
Q

0

CP (u)du + K

−1

A.

0

Integration of Eq. (2.14) is,
Z t
Ct (u) du = U

T

0

Z t
A(u) du.

(2.17)

0Now, substituting Eq. (2.16) into Eq. (2.17) gives,
Z t
Z t
T
−1
Ct (u) du = −U K Q
CP (u)du + U T K −1 A.
0Integration of Eq. (2.13) is,
Z t
Z t
Z t
CP (u) du.
CT (u) du =Ct (u) du + VP
0

工学院大学

(2.18)

00(2.19)

015

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

Substituting Eq. (2.18) into Eq. (2.19) gives,
Z t

T

CT (u) du = −U K

−1

Z t
Q

0

T

CP (u)du + U K−1Z t
CP (u) du

A + VP00= (−U T K −1 Q + VP )Z t

(2.20)CP (u)du + U T K −1 A.

0Dividing Eq. (2.20) by CT (t) results in,Rt

0 CT (u) duCT (t)T= (−U K

−1Rt
Q + VP ) 0

CP (u)du U T K −1 A+
CT (t)CT (t)(2.21)Rt

CP (u)du=α 0+ β.
CT (t)
Equation (2.21) is the LGA. The LGA is such that after some time t∗ , the term,
U T K −1 A
U T K −1 A=,
CT (t)Ct (t) + VP CP (t)(2.22)

becomes constant in time [2, 31]. In turn, Eq. (2.21) becomes a linear relationship between theterms,
Rt

0 CP (u)duCT (t),(2.23).

(2.24)and,Rt

0 CT (u) duCT (t)The slope of the LGA linear relationship is α = −U T K −1 Q + VP , and the y-intercept is
U T K −1 Aβ=
. A linear regression method can then be applied to estimate these regresCT (t)sion parameters, slope and y-intercept.

Equation (2.21) can therefore be expressed simply as,
Y = αX + β,(2.25)

which is a simple linear relationship between X and Y , with a slope of α, and y-intercept β,
where,Rt
X=

0 CP (u)duCT (t)(2.26)
Rt

Y =

工学院大学

0 CT (u) duCT (t)

16,

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

and,
α = −U T K −1 Q + VP
(2.27)U T K −1 Aβ=.
CT (t)The slope of the LGA described in Eqs. (2.21–2.27) is termed DV . DV defines the capacity of tissue receptors to bind the radiotracer [5]. DV is also interpreted as the equilibrium
ratio of the radioactivity concentration in ROI, CT (t) [kBq · cm−3 ], to that in plasma, CP (t)[kBq · mL−1 ]. Equilibrium state is achieved when there is no more irreversible uptake of theradiotracer between the compartments. At this point, all radiotracer flow is due to reversible
uptakes, and all tissue compartments and the arterial plasma compartment are at equilibriumwith each other. In equilibrium state, the radiotracer concentration in the compartments becomes constant, meaning that the flow of the radiotracer into and out of a given compartment
is the same. In other words, equilibrium state is achieved after time t∗ when the y-intercept
term of the LGA has become constant.

The slope and y-intercept term are functions of the kinetic parameters, and it is demonstrated
for the different models as follows:
One-tissue model: α (α1 ) and β (β1 ) are,α1 = −U T K −1 Q + VP
= (−1)(−k2 )−1 K1 + VP=

(2.28)K1+ VP ,
k2

and,

β1 =

=

U T K −1 A
(1)(−k2 )−1 Ct (t)=CT (t)CT (t)Ct (t)+ VP−k2 CT (t)Ct (t)=−k2 (Ct (t) + VP CP (t))
=

工学院大学

1


(t)
−k2 1 + VP CCPt (t)17(2.29).

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

In this case, equilibrium of the LGA is reached when the term CP (t)/Ct (t) in Eq. (2.29) isconstant against time. VP is assumed to be negligible such that the slope and y-intercept termin the expressions in Eqs. (2.28) and (2.29), respectively, are functions of the kinetic parameters
only. Thus, Eq. (2.28) becomes,
α1 =

K1
,
k2

(2.30)

which denotes the DV obtained from a one-tissue compartmental model. Accordingly, Eq.
(2.29) becomes,
β1 = −

1.
k2(2.31)

Two-tissues model: α (α2 ) and β (β2 ) are,α2 = −U T K −1 Q + VP




h
i
−k4
−k4
K
1
  1  + VP
=−
1 1 
k2 k4
−k3 −(k2 + k3 )
0

K1
=
k2



k3
1+
k4

(2.32)

+ VP ,and,β2 =U T U −1 A
CT (t)
− k21k4

h

=

−1
=
k2 k4

工学院大学

1 1

i

−k4



−k4



CF (t)

−k3 −(k2 + k3 )CS (t)(2.33)CT (t)Ct (t)CS (t)(k3 + k4 )+ k2CT (t)CT (t)18

.

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

Three-tissues model: α (α3 ) and β (β3 ) are,α3 = −U T K −1 Q + VP






−k4 k6k4 k6−k4 k6
  K1 
i



1 h
  0  + VP=
1 1 1 k3 k6 k2 k6 + k3 k6
k3 k6




k2 k4 k6
−k4 k5k4 k5−k2 k4 − k4 k5
0

K1
K1
+
=−
k2
k2 k4



k4 k5
k3 −
k6

(2.34)

+ VP ,and,β3 =U T K −1 A
CT (t)




−k
k
k
k
−k
k
C
(t)
4 6
4 6
4 6
 F

h
i



1
  C (t) − k2 k4 k6 1 1 1 k3 k6
 k3 k6 k2 k6 + k3 k6
 S




−k4 k5k4 k5−k2 k4 − k4 k5
CN (t)==

2.2(2.35)CT (t)
1−k3 Ct (t)+k2 k4 CT (t) k21 CN (t)1 CS (t)k5
CF (t) − CS (t) + CN (t)
−
−
.
1+
k6
CT (t)
k6 CT (t)k4 CT (t)Logan graphical analysis using a reference region

Distribution volume ratio and non-displaceable binding potential
Arterial blood sampling to obtain plasma time-activity curve (pTAC), Cp (t), is invasive anduncomfortable for patients [21, 32]. It can also be expensive, which is also undesirable. An
alternative approach is to use a time-activity curve (TAC) of a reference region [3] instead ofpTAC to avoid blood sampling. A reference region has ignorable specific binding sites for theradiotracer and a density of nonspecific binding sites equal to that in the ROI [3]. For the 11 CPiB PET data used in this study, the cerebellum gray matter region was used as the reference
region [33]. This approach of using a reference region TAC in place of pTAC is referred to as

工学院大学

19

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

'

$

&

%

Figure 2.3: One-tissue reference region compartmental model.

reference LGA, and the output parameter is the distribution volume ratio (DVR). DVR is ameasure of the ratio of DV in a receptor-rich region (ROI) to a non-receptor (reference) region,and is used for assessing receptor availability [3].
Z t

Z t

ReplacingCT (u) du and CT (t) in Eq. (2.21) by0Z

and making

CREF (u) du and CREF (t), respectively,
0t

CP (u) du the subject of the formula we get,
0

Z t

1
CP (u) du =
DVREF
0

Z t


CREF (u) du − βREF CREF (t) ,(2.36)

0where the superscript REF is for reference region.

工学院大学

20

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

Z t
Substituting for

CP (u) du in Eq. (2.21) as per Eq. (2.36), we get that,
0

Rt

DV
0 CT (u) du=CT (t)DVREF

Rt

0 CREF (u) du − βREF CREF (t)CT (t)!
+ β.

(2.37)For a two-compartment (one-tissue) reference region shown in Fig. 2.3, we write,
1
βREF = − REF ,k2(2.38)

as per Eq. (2.31). Substituting for βREF in Eq. (2.37) gives,
!
Rt
Rt
REF
C
(u)
du
C
(u)
du
+
C
(t)/k
T
REF
REF
2
0
0
+ β,
= DVR
CT (t)CT (t)
where DVR =(2.39)DV.
DVREF

Replacing k2REF in Eq. (2.39) with an average value, k̄2REF , over all voxels in the reference
region, we get that,Rt
0 CT (u) du= DVRCT (t)Rt

REF
0 CREF (u) du + CREF (t)/k̄2
CT (t)!
+ β,(2.40)

Equation (2.40) is the reference LGA, which uses a reference region TAC instead of pTAC.

Figure 2.4 shows an example of a simulated reference time-activity curve (CREF ) and a ROI
time-activity curve (CT ) in 2.4(a), and a corresponding Logan plot in 2.4(b). The tTACs in Fig.
2.4(a) are averages of 1024 simulated noisy tTACs, which the details will be given in Section
4.1.1. The data plot in the Logan plot in Fig. 2.4(b) achieves a linear relationship at a latertime after the radiotracer administration. These data points, indicated with filled circles, arethe ones used to estimate the Logan slope (DVR in Eq. (2.40)). Time t∗ is then the time point
before the first filled point, which correspond to 30 minutes in Fig. 2.4(a).

DVR can also be interpreted in terms of BPN D , which is related to DVR by,
BPN D = DVR − 1.

(2.41)BPN D compares the concentration of the administered radiotracer in receptor-rich to receptorfree regions [34]. BPN D is proportional to the density of binding sites [3, 34], which denotes theconcentration of the target receptors in ROI. Imaging BPN D is therefore a direct quantification
工学院大学

21

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA'$

Activity [Bq/mL]

104
8

70

7

60

6

50

5
40
4
30
3
20

2

10

1

0

0
0

10

20

30

40

50

60

0

70

5

10

15

20

25

30

35

40

Time [min]

&

(a)

(b)

%

Figure 2.4: An example of a reference time-activity curve (CREF ) and a ROI time-activity curve(CT ) in (a), and a corresponding Logan plot in (b).

of the receptors of interest. BPN D can also be comprehended in terms of its in vitro equivalence0quantity, Bmax KD, which assesses the capacity of the ROI for radiotracer-receptor binding [35].

This quantity is a product of the total number of receptors (Bmax , i.e., free receptor concentra0
tion in the ROI) and their affinity (KD) for the radiotracer [35, 36].

Here, let us take a deeper look into previous methods for bias reduction: Despite
the efficiency and ease of use presented by LGA, further investigations [4, 5, 6, 7] showed thatthe estimates of the LGA slope are negatively biased due to the inherent noise in PET data(tTACs), and that the bias increases with both the noise level and the LGA slope. Due to thesefindings, a number of studies have been carried out to try to reduce this bias. The noise in the
LGA variables is specifically due to the noisy term CT (t) in the denominator of both the LGA
variables, as it can be seen in Eqs. (2.26) and (2.40), presenting a correlated errors-in-variables
problem. The bias is then due to that the OLS regression only accounts for the errors in the
response variable, assuming the predictor variable is error-free. This has been analysed in [4].

One approach [11] for addressing the bias problem is that, since OLS does not account for errors
in both variables, ODR, which minimizes the perpendicular distances from data points to theestimated regression line, does account for errors in both variables and thus would reduce theinfluence of the noise in the estimates of the LGA slope. This approach only removed a part
工学院大学

22

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGAof the bias, and it came at the cost of precision [16]. The reason to not removing much of thebias could be due to that ODR is based on the assumption that the variances of the error invariables are equal, which is generally not the case in practical situations, and specifically not
the case for LGA.

As it is the tradition, assuming that the integral terms in the numerators of Eq. (2.21) arenoise-free, another approach [37] proposed a different method by re-arranging the original LGA
by replacing the noisy term CT (t) in the denominators of Eq. (2.21) by CP (t) and continued asfollows:
Rt

0 CT (u) duCP (t)Rt

CP (u)du U T K −1 A
= (−U K Q + VP ) 0+
CP (t)CP (t)Rt
CP (u)du=α 0+ β.
CP (t)
T−1(2.42)

Multiplying through by CP (t) gives,
Z t

Z t
CT (u) du = α

CP (u)du + βCP (t).

(2.43)00

OLS is then used to estimate DV (slope α) from Eq. (2.43), thus supposedly accounting for
all the noise since the right-hand-side (RHS) of Eq. (2.43) is supposed to be noise-free in thiscase. This is because the unintegrated noisy term, CT (t), no longer appears on the RHS of
Eq. (2.43); it does not appear on the left hand side either. This method did not really bringup about a great improvement, it was shown to be consistent with the classic LGA [37]. ThisZ t
Z t
could be due to that, though these integrals,
CT (u) du and
CP (u) du are assumed to be00

noise-free, they are actually not; they are only less noisy compared to CT (t), but not noise-free.
It was also pointed out in another study [5] that with this formulation, the true equilibrium
state condition must apply before linearity is achieved, and this generally requires a longer timecompared to the classic LGA formulation, with CT (t) in the denominator.

Another approach, multilinear analysis 1 (MA1) [9, 10], is also a mathematical rearrangementof the classic LGA equation into a format that would otherwise minimize the error in the LGA
variables. In MA1, Eq. (2.21) is rearranged to bring the noisy term, CT (t), on one side as,αCT (t) = −
β

Z t

1
CP (u) du +β0Z tCT (u) du.

(2.44)0

OLS can then be used to estimate the regression parameters of Eq. (2.44), and DV would be
工学院大学

23

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA−α1and , as,
β
β
 

β
α
= α.
DV = − −β1found from the ratio of the regression parameters,(2.45)

MA1 reduced the bias but did so at the expense of precison [9, 2]. An extension for MA1
is multilinear reference tissue model 2 (MRTM2). MRTM2 can be obtained by multiplying
through Eq. (2.40) by CT (t) to obtain,
Z t

Z t
CT (u) du = DVR

0

CREF (t)CREF (u) du + 0REF
k̄2
0



+ β ∗ CT (t),(2.46)and making CT (t) the subject of the formula gives,
DVR
CT (t) = − ∗βZ tCREF (t)CREF (u) du + 0REF
k̄2
0



1
+ ∗
β

OLS can then be used to estimate the regression parameters, −Z t
CT (u) du.

(2.47)0DVR
1and ∗ , and DVR can beβ∗βfound as,

 ∗
DVR
β
DVR = − − ∗= DVR.
β1(2.48)The GLLS approach [8] is a temporal smoothing method that attempts to denoise tTACs. In
this approach, each noisy tTAC is separated into two, which are then separately fit to twocompartments (one-tissue) model using GLLS [38]. This method performs well at low noiselevel, but large bias and poor precision has been observed at high noise level [39, 8].

Similar to GLLS, the method of PCA denoising [16] is also a temporal smoothing method inwhich a PCA-based lower dimension linear model is used to denoise tTACs. PCA would reduce
the noise in tTACs by transforming the tTACs into a feature space with the same dimensionality as the original tTACs such that the tTACs could be represented by a reduced number ofdominant features only, while retaining most of the important inherent characteristics of theoriginal tTACs. In this case, all the noise component in the less important PCs would be lost,hence reducing the noise in the reduced dimension data set comprised of the most importantPCs [16]. So the source of the correlated errors in the LGA variables is thus minimised. Thismethod performs considerably well, especially in terms of variance reduction, but there still
remains rooms of improvement in terms bias reduction.

工学院大学

24

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

It is therefore that; OLS is highly biased with moderate precision, ODR removes only a portionof the bias which comes with poor precision. The GLLS, MA1 and MRTM2 methods are significantly effective in removing the bias, but this happens at the expense of precision. PCA quite
maintains precision and removes the bias significantly, but leaving a space for further researches
over a wide range of radiotracers, as well as in terms of determining the appropriate number ofPCs. Not withstanding this tremendous effort to reduce the effect of the noise in PET data on
the estimates of the LGA slope, and though some of the methods significantly reduce the bias,
the bias persists and further studies are necessary.

This study takes a two-step approach to address the bias problem. The first step is to reduce
the noise in the LGA estimates from undenoised tTACs. In the second step, the approachin the first step is employed on denoised tTACs. Specifically, the first step is based on thedemonstrations in [11, 4]. In these studies it can be seen that the regression method used toestimate the Logan slope influences the resulting bias. On this basis, the ”first step” in thisstudy approached the LGA bias problem from the EIV perspective by employing an alternative
linear regression method referred to as LSC to estimate the LGA slope. Unlike OLS or ODR,
LSC accounts for the errors in both the predictor and response variables by minimizing theweighted squared residuals in both the predictor and response variables. Even more, LSC
incorporates the correlation of errors in these variables. Thus, LSC is expected to be appropriate
for correlated errors in the LGA variables. LSC is mathematically superior to both OLS andODR, since it is a general solution to the linear least-squares problem, and also because allother commonly known regression methods, including OLS and ODR, can be derived from LSC
as special cases, under certain assumptions about the statistical properties of the data being
studied. Furthermore, LSC has been shown to improve regression parameter estimates in thefields of geophysics [12] and biogeosciences [40]. LSC was assessed against MRTM2, which is
the most accepted method, and the conventional OLS-based LGA. In the second step, LSC is
combined with tTACs denoising techniques, PCA and CorrCA.

Ordinary least-squares
Since LGA is a linear equation. Let us first refresh on OLS and ODR regression methods.
OLS regression is a statistical method used to estimate the relationship between one or morepredictor variables and a response variable, by estimating the unknown parameters in a lin工学院大学

25

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

'

$

&

%

Figure 2.5: A schematic figure for OLS.

ear regression model. OLS estimates this relationship by minimizing the total squared errors[41, 42, 43], i.e., the sum of the squared differences between the observed and the predictedvalues of the response variable. Geometrically, this is the sum of the squared distances parallelto the axis of the response variable, between each scattered data point and the corresponding
point on the estimated regression line. This interpretation can be seen in Fig. 2.5. Since the
LGA is a bivariate system, i.e., one predictor variable and one response variable, henceforward
OLS is presented based on a bivariate system.

Let us define X = {Xi | i = 1, 2, ..., n} and Y = {Yi | i = 1, 2, ..., n} as the predictor and response
variables, respectively. With reference to Fig. 2.5, the OLS line of best fit, y = a + bx, minimizes
the quantity,
D=

n
X

d2i

i=1

(2.49)nX=(Yi − yi )2 ,
i=1

工学院大学

26

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

subject to,yi = a + bXi

fori = 1, 2, · · · , n,(2.50)

where {yi | i = 1, 2, ..., n} are the estimated data points, lying on the estimated straight line.
Substituting Eq. (2.50) into Eq. (2.49) gives,D=nX(Yi − bXi − a)2 .

(2.51)

i=1

We therefore wish to solve for b and a by minimizing Eq. (2.51).

To minimize Eq. (2.51), consider that,D=nX(Yi − bXi − a)2 = D(b, a).

(2.52)

i=1

We then differentiate Eq. (2.52) with respect to b and a to get,
n

X
∂d2
= −2
(Yi − bXi − a) Xi ,∂b

(2.53)

i=1

and
n

X
∂d2
= −2
(Yi − bXi − a) .
∂a

(2.54)

i=1

Setting Eq. (2.54) to zero we get,nX(Yi − bXi − a) = 0

i=1(2.55)=⇒ a = Ȳ − bX̄,
where X̄ =

n
X
i=1

Xi and Ȳ =

n
X

Yi .

i=1

Setting Eq. (2.53) to zero we get,nX(Yi − bXi − a) = 0, and substituting for a as per Eq.

i=1(2.55) gives,

工学院大学

27

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGAnXYi − bXi − Ȳ + bX̄ = 0

i=1

=⇒ b

n
X

n
 X

Xi − X̄ =Yi − Ȳ

i=1

i=1(2.56)

P
Yi − Ȳ=⇒ b = P
Xi − X̄

P
∴

b=



Yi − Ȳ Xi − X̄.
2P
Xi − X̄

Thus, Eqs. (2.56) and (2.55) respectively give the OLS estimates of slope b and y-intercept aof the straight line in Eq. (2.50).

We can see in Eq. (2.49) that OLS only minimizes the squared residuals of the response variable. This is based on the assumption that only the response variable is contaminated withmeasurement errors. However in practice, this is not always the case, and specifically for LGA
as described earlier, both variables are contaminated and the contamination in both variables
need to be accounted for.

In order to address the issues encountered with OLS, a variety of regression methods referred
to as EIV have been developed. Before discussing the LSC method which is employed in thisstudy, let us have a look at one of the EIV methods, ODR, which has been previously applied
to LGA [11].

Orthogonal distance regression
ODR (also referred to as PCA linear regression in some literature, eg., [6]) is an EIV regression
method that was developed to account for errors in both the predictor and response variables byminimising the perpendicular distances from the estimated regression line to the observed data
points [44, 45], unlike OLS which minimizes the vertical distances. Fig. 2.6 shows a schematic
diagram for ODR.

In ODR, both Xi and Yi measurement are considered to be noisy. Therefore, ODR estimates
工学院大学

28

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

'

$

&

%

Figure 2.6: A schematic figure for ODR.
the best fit line, y = a + bx, by minimizing the residuals in both the Xi and Yi measurements
as represented in the equation,
D=N
X(Xi − xi )2 + k(Yi − yi )2 ,(2.57)

i=1

subject to,yi = bxi + a,

(2.58)

where {(xi , yi ) | i = 1, 2, ..., n} are the estimated points sitting on the estimated straight line, and22
k is the ratio of variances of the errors in the measured data given by σxi/σyi. Note here thatEq. (2.58) differs from Eq. (2.50), in that Eq. (2.58) has xi instead of Xi which is in Eq. (2.50).

Substituting Eq. (2.58) in Eq. (2.57) gives,
D=NX(Xi − xi )2 + k(Yi − bxi − a)2 .

(2.59)

i=1

The respective solutions of the slope and y-intercept from Eq. (2.59) are,q
2
(kSyy − Sxx ) + (Sxx − kSyy )2 + 4kSxy
b=
,
2kSxy
工学院大学

29

(2.60)KOGAKUIN UNIVERSITYPK SHIGWEDHA

PET, COMPARTMENTAL MODEL AND LGA

anda = Ȳ − α0 X̄,
where Sxx =X(Xi − X̄)2 , Syy =X(Yi − Ȳ )2 and Sxy =(2.61)
X(Xi − X̄)(Yi − Ȳ ).

In Eq. (2.60), if we assume that σxi = σyi , then k = 1, and Eq. (2.60) becomes,q2(Syy − Sxx ) + (Sxx − Syy )2 + 4Sxy
.
b=
2Sxy

(2.62)

Equations (2.62) and (2.61) are the ODR solutions of the slope and y-intercept, respectively.
Thus, ODR is mainly based on this assumption that the variances of the errors in the variables
are equal. Clearly this is not always the case in practical situations, that is why ODR regression
does not make the best regression method. Next, I present LSC linear regression method
according to [13, 14, 12, 46] and show how it can be used to estimate the LGA slope.

工学院大学

30

KOGAKUIN UNIVERSITY

Chapter 3

Reducing the bias in the BPN D
estimates using least-squares cubic
linear regression
3.1

Least-squares cubic linear regression

LSC is an EIV-based regression method. Building upon the work of Demming [47], LSC wasdeveloped by Derek York in (1966 [13]), in which he considered the errors in both the predictor
and response variables but not the correlation of the errors in variables. He later (1969 [14])extended it to take into account the correlation of measurement errors in variables. Further,in 2004, York and colleagues [12] followed up, and presented compact equations for the slopeof LSC, and unified the the expressions of the standard errors of the slope and y-intercept
with maximum likelihood expressions. LSC fully considers the errors in both the predictor andresponse variables [12]. This is achieved by minimizing the sum of weighted squared errors inboth variables, and by including the correlation of these errors [12]. Under specific assumptions,
LSC can be shown to reduce to either OLS or other common regression methods (see Appendix
A), rendering LSC to be a general solution to the linear regression problem [13, 14].

To discuss the mathematical context of LSC, consider the observed predictor and response
variables, X = {Xi | i = 1, 2, ..., n} and Y = {Yi | i = 1, 2, ..., n}, respectively. The LSC slope is

31

PK SHIGWEDHA

LSC

obtained by minimizing the equation below [12];
S=

n
X


w(Xi )(xi − Xi )2 − 2ri αi (xi − Xi )(yi − Yi ) + w(Yi )(yi − Yi )2

i=1

1
,
1 − ri2(3.1)

where (xi , yi ) are the estimated points of (Xi , Yi ) sitting on the estimated straight line,pand αi = w(Xi )w(Yi ). Equation (3.1) includes the residuals in both variables [(xi − Xi ),(yi − Yi )], the weights [w(Xi ), w(Xi )], and the correlation of the errors in variables (ri ).

Deriving Eq. (3.1)
The derivation of Eq. (3.1) is demonstrated in [46] as follows:
For the respective ith observed and estimated points (Xi , Yi ) and (xi , yi ), let us denote
the residuals column vector by,
viT = [u(Xi ), v(Yi )] ,(3.2)

where, superscript T denotes a transpose, andu(Xi ) = xi − Xi
(3.3)
v(Yi ) = yi − Yi .
We can represent the column vector of residual vectors as,vT = [v(X1 ) v(Y1 ) | v(X2 ) v(Y2 ) | · · · · · · | v(Xn ) v(Yn )]
(3.4)


= v1T , v2T , · · · vnT .

Let s2 (Xi ), s2 (Yi ) and s(Xi Yi ), respectively, denote the variances and covariances
associated with the observed pair (Xi , Yi ). The block-diagonal variance matrix is defined

工学院大学

32

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

as,



s2 (X1 )s(X1 Y1 )
0
0




 s(X1 Y1 ) s2 (Y1 )
0
0




0
0
s2 (X2 ) s(X2 Y2 )




0
0
s(X2 Y2 ) s2 (Y2 )


Q=
..
..
..
..


.
.
.
.



..
..
..
..

.
.
.
.




0
0
0
0



0
0
0
0



···

···

···

···

···

···

···

···

..

.
..

.

···

···

···

···

0

0






0
0




0
0




0
0



..
..


.
.



..
..

.
.



s2 (Xn ) s(Xn Yn ) 2s(Xn Yn ) s (Yn ) Q1 0 · · · 0 Q2 · · ·

=
 .
..
..
 ..
.
.



0
0 ···

0 



0 

,
.. 
. 



Qn
(3.5)

in which the 2-by-2 sub-matrices, Q1 , Q2 , · · · , Qn , are symmetric.

If we refer to

X1 , X2 , · · · XN as a finite population of N quantities, the population mean and variance are respectively given by,
N

1 X
Xi
µ=
N
i=1

N

1 X
and σ =(Xi − µ)2 .
N
2

(3.6)

i=1

Given a sample of n quantities from an infinite population, the unbiased estimates of themean and variance are respectively,
n

X̄ =

1X
Xi
n

n

and s2 (X) =

i=1

2
1 X
Xi − X̄ .
n−1

(3.7)

i=1

The standard deviation is given by the positive square-root of the variance, i.e.,s2nX
Xi − X̄
s(X) =
.
n−1

(3.8)

i=1

For a pair of two samples (Xi , Yi ) both of size n, an unbiased estimate of covariance is,ns(XY ) =
1 X
Xi − X̄ Yi − Ȳ .
n−1

(3.9)

i=1

工学院大学

33

KOGAKUIN UNIVERSITYPK SHIGWEDHALSC

The block-diagonal weight matrix of the measurements denoted by W is defined by,


0
···
 W1



 0
W2 · · ·W = Q−1 
 .
..
..
 ..
.
.



0
0
···

0
0
..
.
Wn







,






(3.10)−1−1
in which the 2-by-2 sub-matrices, W1 = Q−1
1 , W2 = Q2 , · · · , Wn = Qn , aresymmetric. This means that the weight of measurements are given by the inverse of thevariances.

From Eq. (3.5) we can extract the ith variance matrix Qi as,
Qi = s2 (Xi )s(Xi Yi ) ,
s(Xi Yi ) s2 (Yi )(3.11)and the inverse of the ith variance matrix will give the ith weight matrix Wi as,
Wi = Q−1
=
i


1

2
2
2
s (Xi )s (Yi ) − s (Xi Yi ) 
s2 (Yi )−s(Xi Yi ) .

−s(Xi Yi )s2 (Xi )(3.12)

The correlation (ri ) between a pair of measurements (Xi , Yi ) is given by,ri =s(Xi Yi ),s(Xi )s(Yi )(3.13)

which gives,1s2 (Xi )s2 (Yi ) − s2 (Xi Yi )

工学院大学

=

34

1
1 − ri2s2 (Xi )s2 (Yi ).

(3.14)KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Substituting Eq. (3.14) into Eq. (3.12) gives,






1


s2 (Xi )1 Wi =1 − ri2 



 −
ri

i
− s(Xir)s(Yi)

1s2 (Yi )s(Xi )s(Yi )






.






(3.15)

Expressing the weights as the inverses of the unbiased estimates of the variances we get,w(Xi ) =1
s2 (Xi )and w(Yi ) =1
s2 (Yi ).

(3.16)

We can then write Eq. (3.15) as,






 w(Xi ) −ri (αi ) 

1 
,

Wi =


2
1 − ri −ri αi w(Yi )(3.17)

1.
s(Xi )s(Yi )(3.18)

where,αi =pw(Xi )w(Yi ) =The least-squares function S in Eq. (3.1) is equal to the sum of weighted squares ofresiduals, and can be defined as a matrix product,
S = vT Wv.

(3.19)By Eqs. (3.4) and (3.10), we can write S as below,
S = v1 T W1 v1 + v2 T W2 v2 + v3 T W3 v3 + · · · + vn T Wn vn ,

工学院大学

35

(3.20)KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

and using Eq. (3.17) we can write,









 w(Xi ) −ri (αi )   v(Xi ) 




1


v i T Wi v i =[v(Xi ) v(Yi )] 




2
1 − ri






v(Yi )−ri αi w(Yi )

=

=

1 
w(Xi )(v(Xi ))2 − 2ri αv(Xi )v(Yi ) + w(Yi )(v(Yi ))2
1 − ri2
1 
w(Xi )(xi − Xi )2 − 2ri α(xi − Xi )(yi − Yi ) + w(Yi )(yi − Yi )2 .
1 − ri2(3.21)

Summing over i = 1, 2, · · · , n gives,
S=

n
X

1 
2
2
2 w(Xi )(xi − Xi ) − 2ri α(xi − Xi )(yi − Yi ) + w(Yi )(yi − Yi )
1
−
r
i
i=1
(3.22)

as in Eq. (3.1).

Solving for the least-squares cubic line of best fit
The line of best fit y = a + bx minimizes the sum of the weighted squares of residuals [13, 14].
The adjusted/estimated points (xi , yi ) of the measurements (Xi , Yi ) are required to lie on theestimated straight line, and therefore the regression parameters a,

and b can be obtained byminimizing S (Eqs. (3.22) and (3.1)) subject to the constraints that,yi = a + bxi

fori = 1, 2, · · · , n.

(3.23)

It is shown in [46] how this task can be achieved using the Lagrange multipliers optimization.
The function to be minimized is the Lagrangian L written as,
L=S+2

X

λi (a + bxi − yi ) ,(3.24)iwhere λi are unknown Lagrange multipliers. The Lagrangian, L, has five variables, λi , a,b, xk , and yk , and can be minimized by setting its first partial derivatives to zero as below,
∂L
∂L
∂L
∂L
∂L
=
=
=
=
= 0.
∂λi
∂a
∂b
∂xi
∂yi

工学院大学

36

(3.25)KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Working out the zero-equated partial derivatives in Eq. (3.25) gives,X∂L=2(a + bxi − yi ) = 0
∂λi
i

X
∂L
=2
λi = 0
∂a
i

X
∂L
=2
λi x = 0∂b(3.26)

i

X 1
X
∂L
=
{2w(X
)(x
−
X
)
−
2r
α
(y
−
Y
)}
+
2
λi b = 0
i
i
i
i
i
i
i
∂xi
1 − ri2
i

i

X
X 1
∂L
{2w(Y
)(y
−
Y
)
−
2r
α
(x
−
X
)}
−
2
λi = 0.
=
i
i
i
i
i
i
i
∂yi
1 − ri2ii

In order, the expressions in Eq. (3.26) gives,
yi = a + bxi(3.27)

λi = 0(3.28)

λi x = 0(3.29)w(Xi )(xi − Xi ) − ri αi (yi − Yi ) + λi (1 − ri2 )b = 0(3.30)

w(Yi )(yi − Yi ) − ri αi (xi − Xi ) − λi (1 − ri2 ) = 0,

(3.31)

for each i.

Respectively solving for xi − Xi and yi − Yi from Eqs. (3.30) and (3.31) gives,
xi − Xi =1
1(ri αi λi − bw(Yi )λi ) =⇒ xi = Xi + 2 (ri αi λi − bw(Yi )λi )2αi
αi(3.32)

yi − Yi =1
1(w(Xi )λi − bri αi λi ) =⇒ yi = Yi + 2 (w(Xi )λi − bri αi λi ) .
αi2
αi(3.33)

Substituting Eqs. (3.32) and (3.33) into Eq. (3.27) and solving for the Lagrange multipliers
gives,
λi = Wi (a + bXi − Yi ),(3.34)

where Wi is the weight function given by,
Wi =
工学院大学

αk2
.
b2 w(Yi ) + w(Xi ) − 2bri αi37(3.35)KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Substituting Eq. (3.34) into Eqs. (3.28) and (3.29), respectively, gives,X

Wi (a + bXi − Yi ) = 0(3.36)

xi Wi (a + bXi − Yi ) = 0.

(3.37)

andX

Expanding Eq. (3.36) and solving for a gives,

∴

P
P
Wi Xi
Wi Yi
−b P
a= P
Wi
Wi

(3.38)

a = Ȳ − bX̄,

(3.39)

where (X̄, Ȳ ) is the Centroid of the data:P
Wi Xi
X̄ = P
,
Wi

P
Wi Yi
Ȳ = P.
Wi(3.40)

Equations (3.38) and (3.39) mean that, according to the LSC approach, the best straight line
goes through the centroid/center of gravity of the data set given by (X̄, Ȳ ) in Eq. (3.40).

Let us define the centroidal coordinates (Ui , Vi ) as,
Ui = Xi − X̄,Vi = Yi − Ȳ .

(3.41)

Then we can write Eq. (3.39) asa + bXi − Yi = a + b(Ui + X̄) − (Vi + Ȳ )(3.42)= a + bUi + bX̄ − Vi − Ȳ
= bUi − Vi
Respectively substituting Eq. (3.42) into Eqs. (3.34) and (3.37), we get,λi = Wi (bUi − Vi ) = 0(3.43)

X(3.44)

and工学院大学

Wi (bUi − Vi ) = 0.

38

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Substituting Eqs. (3.32), (3.42), (3.43), and (3.44) into Eq. (3.37), gives,


X
X W2
2
ri
2
2
2
3
i
Ui +
U −b
Wi
V i Uibw(Xi ) i
αiw(Xi )(3.45)
−b

X




 X
Wi
Wi ri
2 2Wi ri2
W i UkUi V i −Vi = 0.
V+
W i Ui −αiw(Xi ) i
αi

Equation (3.45) is referred to as the Generalized least-squares cubic regression [14, 12].
Being a general solution to the least-squares problem means that all other known regression
methods, like ODR and OLS can be obtained from Eq. (3.45) as special cases, under specificassumptions. In addition, the solution obtained by York in [13] in which it is assumed thatthe errors in the variables are uncorrelated can also be obtained from Eq. (3.45) by setting ri
to zero. The “cubic” in “Generalized least-squares cubic regression” refers to that Eq.
(3.45) is a cubic (or rather pseudo-cubic, as explained below) equation of b.

The pseudo-cubic equation:
Equation (3.45) is not really a cubic equation of b because the qantity Wi (Eq. (3.35))
and the centroidal coordinates (Ui , Vi ) (Eq. (3.41)) are functions of b. However, thispseudo-cubic equation can be reduced to a truc cubic equation by substituting in an
approximate value for b in the weight function Wi . From there, it can then be solved
directly as a cubic equation, giving three roots for b.

A concise solution for b is given in [12]. It was obtained by first reducing Eq. (3.45) to aquadratic equation as,
2

b

X

Wi2



Ui V irU 2− iw(Xi )
αi


+b

X

Wi2



Ui2
Vi2
−
w(Yi ) w(Xi )


−

X

Wi2



Ui V irV 2− iw(Yi )αi= 0,(3.46)and subsequently to a linear equation as,
b

X

Wi2



Ui
bVi
bri Ui+−w(Yi ) w(Xi )
αi


Ui −

X

Wi2



Ui
bVi
ri Vi
+
−
w(Yi ) w(Xi )αi

The slope coefficient, b, can then be obtained from Eq. (3.47) as,
P 2  UibViri Vi
Wi w(Yi ) + w(X
−
Vi
αi
i)

 .
b= P
Ui
bVi
Wi2 w(Y+ w(X− brαi Ui i Ui
i)
i)

工学院大学

39


Vi = 0.

(3.47)(3.48)KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Further, as shown in [46], Eq. (3.47) can be written as,




X
X
Ui
ri Vi
bVi
bri Uib
W i Ui W i−b
W i Ui W i+−w(Yi ) w(Xi )αiαi(3.49)=X
Wi Vi Wi

Ui
bVi
ri Vi
+
−
w(Yi ) w(Xi )
αi


−b

X


W i Ui W i

ri Vi
αi


,

which can be further simplified to give the expression,

 X
X
Ui
bVi
ri
b
Wi Ui Wi
Wi V i Wi+− (bUi + Vi ) =w(Yi ) w(Xi ) αi

 (3.50)
Ui
bVi
ri Vi
(bUi + Vi ) .
+−w(Yi ) w(Xi )αi

The LSC slope parameter, b, can then be expressed in simplified terms as,P
Wi βi ViPb=.
Wi βi Ui

(3.51)

where βi is given by,
β i = Wi

bVi
ri
Ui
+
− (bUi + Vi )w(Yi ) w(Xi )αi.

(3.52)

LSC solution for the slope, b, and y-intercept, a:
P
Wi Bi Vi
b= P
W i B i Ui

from Eq. (3.51)

a = Ȳ − bX̄

from Eq. (3.39)

Respective variances of a and b:
σb2 = P1
Wi u2i(3.53)

1
σa2 = P+ x̄2 σb2
Wi

Estimated/adjusted points (xi , yi ):xi = X̄ + βi
(3.54)
yi = Ȳ + bβi

工学院大学

40

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Residuals (xi , yi ):v(Xi ) = xi − Xi = βi − Ui(3.55)
v(Yi ) = yi − Yi = bβi − Vix̄ in Eq. (3.53), and the corresponding ȳ are respectively given by,P
Wi xix̄ = PWiP
Wi yi
and ȳ = P,
Wi(3.56)

ui = xi − x̄and vi = yi − ȳ.

(3.57)

where,

Here let us again note that Wi itself (which is a function of b), and βi , Vi and, Ui , which are allfunctions of Wi (and therefore functions of b) appear in Eq. (3.51). Therefore an appropriate
approximate value for b shoul be substituted in Wi so that Eq. (3.51) can be used to estimateb. In this study, the approximate value which was substituted in Wi was estimated by OLS.
For improved estimates of b, the estimates were iterated. A complete mathematical feature of
LSC is provided in [46].

工学院大学

41

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

'

$

&

%

Figure 3.1: An illustration of LSC regression: LSC for a general Y vs X. The observed point(Xi , Yi ) is projected onto the estimated point (xi , yi ), making a residual di .

工学院大学

42

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

'

$

&

%

Figure 3.2: An illustration of LSC regression: LSC visualized in terms of the LGA variables.
The main point is that the direction along the minimization distance is not restricted to being
either vertical(as is for OLS)R
R or perpendicular (as is for ODR), but may vary depending on theCT (u) duCT (u) duvalue of R
ξ/ξ = R, as implied in Eq. (3.58).
CP (u) duCP (u) du

工学院大学

43

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Analysing least-squares cubic regression in terms of the Logangraphical analysis
Figure 3.1 shows a schematic diagram for a general LSC regression of Y vs X, and Fig. 3.2
expresses LSC specifically in terms of the LGA variables.

To explain LSC in specific terms of LGA, let us consider the LGA in Eq. (2.21). If theZZ
two integral terms,
CP (u) du andCT (u) du , are free from measurement errors, then allRR
CT (u) duCp (u) du
the errors in the LGA variables,, and, in Eq. (2.21) are entirely due
CT (t)CT (t)
to CT (t). Accordingly, the errors in these two LGA variables will be perfectly correlated, i.e.,ri = 1. The assumption that the integral terms are noise-free is based on the fact that integrated
measurements has less noise than the discrete data. To simplify
the explanation, we assumeRCP (u) du
this is the case. Given the above, lets us say the value ofat point H in Fig. 3.2
CT (t)R
CT (u) du
has an error of ξ due to CT (t) error, then, according to York in [13], the value ofCT (t)R
CT (u) duξ due to that same CT (t) error. Now, the observed point H iswill have an error of R
CP (u) du
shifted from the true point H ∗ along a line with a slope given by,
∆Y
=
∆X

Rt
CT (u) duR0t
ξC
0 P (u) du

ξ
Rt

= R0t

CT (u) du

0 CP (u) du(3.58).

In a less compact format than Eq. (3.55), the LSC residuals can be expressed as,xi − Xi =Wi (a + bXi − Y ) (ri αi − b w(Yi ))w(Xi )w(Yi )(3.59)

Wi (a + bXi − Y ) (w(Xi ) − bri αi )y i − Yi =.
w(Xi )w(Yi )Using Eq. (3.59), the slope of the straight line connecting points H and H ∗ in Fig. 3.2 is givenby,dslope =工学院大学

Yi − yiXi − xi=y-residual
x-residual=w(Xi ) − bri αi.
ri αi − b w(Yi )44(3.60)KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Under the assumption that the errors are perfectly correlated, we have, r = 1, and Eq. (3.60)
becomes,pw(Xi ) − b w(Xi )w(Yi )dslope = pw(Xi )w(Yi ) − b w(Yi )(3.61)s=w(Xi ).
w(Yi )

From Eq. (3.58) we have that,
Rt

Rt
C (t) duCT (t) duerror in X measurement
∆Y
=
= R0t
,
=
Rt
error in Y measurement
∆X
C
(t)
du
0 CP (t) duPerror inT
error in 0 C(t)(3.62)0C(t)and by using the definition of the weights as the inverse of the squared errors, we get that,


Rt
Rt

CP (t) dus#
Rt
1/ error in 0 C(t)0 CT (t) du
error
in

C
(t)
du
w(Xi )TC(t) = R0,(3.63)
=
Rt
Rt
t
 =
CP (t) duw(Yi )
C
(t)
du
C
(t)
du
0
T
P
0
error
in
0
i
1/ error in
C(t)
C(t)
i
i

i.e.,
s

#
Rt
C
(t)
du
w(Xi )T.
= R0t
w(Yi )CP (t) du0(3.64)iEquation (3.61) serves to inform that the residuals are not necessarily, vertical, horizontal or
perpendicular to the estimated straight line. They can be at any angle which is defined by theirslope. We see in Eq. (3.61) that the slope of the residuals is dependent on the weight functions.
This means that the weights will serve to determine the direction in which the observed measurements should be adjusted.

To put it into perspective, for the unweighted OLS, w(Yi ) = 1 and w(Xi ) = 0, and therefore
dslope = inf. This means that the residuals of the OLS regression are in the vertical direction,
which is exactly the case for OLS, minimizing only the errors in the Y measurements.

Equation (3.64) shows how the weights —when assumed to be the inverse of measurement
errors— relate to the integral terms in the LGA variables. Given that the errors in measured
data are practically unknown, the aim is to define the weight functions such that they satisfy
the relationship in Eq. (3.64).

工学院大学

45

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

The weights [w(Xi ), w(Yi )] and the correlation of the errors ri
The variances in the measurements are not known, and therefore the weights cannot be directly
calculated by,w(Xi ) =1
s2 (Xi )andw(Yi ) =1s2 (Yi ).

(3.65)

This study adopted the method of “relative error regression” [48, 49]. In relative error regression,
X
instead of minimizing the squared-sum of the absolute errors,(yi − Yi ), the squared-sum ofXthe relative errors,yi − Yi )/Yi2 , is minimized. In statistical prediction, relative errors are
considered to be more informative than absolute errors [49]. Expressing the relative errors inthis manner can help adjust measurements in a weighted linear regression with the weights
given by, 1/Yi2 [48, 49]. After adopting this technique for LSC, the X and Y weights can berespectively given as,w(Xi ) =

1
Xi2

and

w(Yi ) =1.
Yi2(3.66)

These weight functions suffices for the LSC-based LGA since they satisfy Eq. (3.64) as below,
s
s
1/Xi2
Yi
w(Xi )==
2
w(Yi )Xi
1/Yi
R
CT (ti ) du/CT (ti )R=CP (ti ) du/CT (ti )(3.67)R
CT (ti ) du=R.
CP (ti ) du

The correlation of the errors were estimated by the correlations between the variables themselves. These correlations were estimated by the correlations of determination. The correlationof determination is a number between 0 and 1, and it determines the proportion of variation in
the response variable that can be explained by the predictor variable [50, 51]. To iterate, theslope obtained in the first run was substituted in Wi and βi . The iteration process was repeateduntil the relative differences between consecutive estimates of the slope, (bj − bj+1 ) /bj+1 , wereless than 10−10 .

In summary, as per [12], here is how one would carry out the estimation of the regression
parameters by LSC;

工学院大学

46

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

1. Begin by obtaining an approximate initial value for b. In this study, it was obtained by
OLS.
2. Approximate the weights, w(Xi ) and w(Yi ), for each point, and the correlation between
the errors in measurements, ri . In this study, the weights were obtained by Eq. (3.66),
and ri was estimated by the correlation between the measurements themselves.
3. Using the weights, the initial approximate for b and the correlation, ri , evaluate Wi foreach point.
4. Calculate X̄ and Ȳ by Eq. (3.40), calculate Ui and Vi by Eq. (3.41), and then calculate
βi by Eq. (3.52), for each point.
5. Using Wi , Ui , Vi , and βi , calculate b by Eq. (3.51).
6. Now use the new b in step 5 as the initial estimate and repeat steps, 3, 4, and 5 until
consecutive estimates of b converges to some desired condition.
7. With this final value of b, together with the final X̄ and Ȳ , calculate a by Eq. (3.39).
8. For each (Xi , Yi ), calculate the adjusted/estimated points, xi and yi , by Eq. (3.54).
9. Use the estimated points, (xi , yi ), and Wi to calculate (x̄, ȳ) by Eq. (3.56), and thencalculate (ui , vi ) by Eq. (3.57).
10. Using Wi , x̄ and ui , calculate σb , and then σa by Eq. (3.53).

3.1.1

Simulation studies

Simulation data design
In order to assess the performance of LSC against that of MRTM2 and conventional LGA
when estimating BPN D , a set of PET data of a two-tissue compartmental model was simulatedusing a clinically measured pTAC. The simulated data mimicked a well-known radiotracer,
11

C-carfentanil (CFN), which is a reversible radiotracer that binds to µ-opioid in the brain.

The kinetic parameters, K1 , k2–4 , and the non-displaceable distribution volume (VN D ) used for
simulations were adopted from a range of values used in a previous study [52];• 11 values of BPN D were set in the range of [0.0, 3.0], covering a range of [0.0, 0.35] for
k3 .
工学院大学

47

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC• [K1 k4 ] = [0.1835 mL cm−1 min−1• [k2 k3 ] = [K1 /VN D• VN D

0.115/min].

BPN D · k4 ].

= 1.59 mL cm−1 .

Simulated data were made for a 90 min, 33-frame scan (6 × 0.1667 min, 3 × 0.5 min, 5 × 1 min,5 × 2.5 min, and 14 × 5 min). The reference region was formed with a zero noise level, zeroBPN D , and delivery and clearance rates equal to those of the target tissues, due to its physiological assumptions. Using these kinetic parameter values, 11 noise-free tTACs corresponding
to the 11 BPN D values were formed. Statistical noise was added to these noise-free tTACs, and1024 noisy tTACs were simulated as a slice of 32-by-32 pixels for each noise-free tTAC. The
BPN D values were then re-estimated from the noisy tTACs by LSC, MRTM2, and conventional
LGA, and the results were compared. The true value for the dissociation rate of the reference
region (k2R ), which was used to produce the simulation data, was used for all methods.

The noise was assumed to have a zero-mean Gaussian distribution with a variance proportional
to the true tTAC; this is the same to the noise model used in previous studies [6, 8, 16, 11]. The
noise scaling factor was set such that the actual magnitude of the noise level in the simulated
tTACs varied within a range of 0 to 30%; this range covered the actual noise observed in voxelbased PET data. The magnitudes of the noise level were calculated by the percentage of theratios of the “standard deviation of the noisy tTACs” to the “mean of the noise-free tTAC” forthe portion of the time range used for the BPN D estimations.
Results
Figure 3.3 shows the simulated noise-free (thick lines) and their corresponding noisy (thin lines)
tTACs for the 11 BPN D values. The 12th blue curve at the bottom is the curve of the reference
region. Out of the the 1024 simulated noisy tTACs, only one is shown in Fig. 3.3 for eachnoise-free tTAC.

The results in Fig. 3.4 compares the percentage bias (difference from 100%) in the BPN D values estimated by three methods, LSC, MRTM2 and OLS, from the noisy tTACs. The standard
deviations (in percentages) are shown as error bars. In these results, LSC estimates are theleast biased along the entire range of BPN D , with up to a maximum of 12% and 28% bias
工学院大学

48

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

'

$

&

%

Figure 3.3: Simulated tTACs. The 12 thick curves denote the noise-free tTACs. Eleven of these
correspond to the 11 BPN D values. The remaining one (lowermost) represents the reference
region. The thin curves denote the noisy tTACs. Only one noisy tTAC is shown out of the 1024
simulated noisy tTACs for each of the noise-free tTACs.

difference in comparison with MRTM2 and conventional OLS-based LGA, respectively. The
error bars show that while the conventional LGA estimates have the smallest variability, LSC
estimates have smaller variability than MRTM2 estimates. With the above, it can be seen that
LSC outperforms MRTM2 in terms of both bias and variance, i.e., LSC estimates have both
smaller bias and variance than MRTM2 estimates.

Figure 3.5 shows the distributions of the errors obtained from LSC- and OLS-based LGAs.
The errors for OLS refers to the y-residuals. The errors for LSC are calculated as a sum of
Eqs. (3.27–3.31), which makes the condition for the minimization term to equate to zero. The
histograms in Fig. 3.5 are for residuals of only 8 random of the simulated tTACs. In this figure,it can be seen that LSC estimates make a better normal distribution —both in terms of the

工学院大学

49

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

'

$

&

%

Figure 3.4: Comparisons of the percentage bias and variability in the estimated BPN D values.
The percentage bias and variability in the BPN D estimated from the noisy tTACs by LSC,
MRTM2, and conventional LGA, are compared for the 11 values of true BPN D .

bell shape and in terms of narrowness— in comparison to those of OLS. Figure 3.6 shows thehistograms of the errors of all 1024 simulated tTACs for BPN D = 3 —a further confirmation of
the results in Fig. 3.5.

3.1.2

Clinical 11 C-PiB PET data studies

Study design
Two 11 C-PiB PET data sets, of which one is Aβ-negative and the other is Aβ-positive, wereused for this analysis. The two participants underwent a 11 C-PiB PET scan for 70 minutes with
25 frames (6 × 0.1667 min, 3 × 0.3333 min, 2 × 1 min, 2 × 3 min, and 12 × 5 min). The scan
was carried out according to the standard ADNI 11 C-PiB PET procedure. The study protocol
was approved by the Ethics Committee of Kindai University Hospital, and written informed
工学院大学

50

KOGAKUIN UNIVERSITY

工学院大学51%$

and the lower panel shows those of OLS.

Figure 3.5: Comparison of the distributions of LSC against OLS errors. The upper panel shows the distributions of the errors estimated by LSC,&'

PK SHIGWEDHALSC

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

consent was obtained from the participants.

The data were collected in arrays of 128 voxels × 128 voxels × 47 slices, with voxel sizes of
2.1 × 2.1 × 3.4 mm3 . BPN D parametric images – reflecting Aβ deposits in the brain regions –were then generated using three algorithms, LSC, MRTM2, and conventional OLS-based LGA,
implemented in MATLAB (The MathWorks, Inc., Natick, MA, United States).

Time t∗ (30 min.), which was used for both synthetic and actual data, was graphically determined from simulated data.

For the actual data, the values of k2R that were used for all methods were estimated by MA1.
MA1 allows estimation of k2R with little bias [10]. Specifically for each slice, the average value,k̄2R , calculated from the | k2R | of all voxels in the reference region, was used.
Results
Figure 3.7 shows the BPN D parametric images obtained by the three methods, LSC, MRTM2,
and conventional OLS-based LGA. These images are generated from 11 C-PiB PET brain images, and therefore they reflect Aβ deposits in the brain regions. The left panel shows imagesof three slices of the Aβ-negative patient, obtained by each of the three methods, whereas theright panel shows the images of the Aβ-positive patient.

In Fig. 3.7, the images obtained by LSC have significantly higher total BPN D estimates thanthose of conventional OLS-based LGA for both patients. The images obtained by MRTM2 havethe highest total BPN D estimates. For MRTM2 images, however, considerably higher levels ofnoise are observed. These results are consistent with those of CFN observed in the simulated
data.

Table 3.1 shows the average computational time taken by each method for the corresponding
sets of three images in Fig. 3.7. LSC took moderate computational time, one order of magnitude
higher than conventional OLS-based LGA, for both patients’ images. This was expected of LSC
because it is an iterative procedure. MRTM2 took the longest computational time amongst thethree methods, two orders of magnitude higher than conventional OLS-based LGA, for both

工学院大学

52

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

Table 3.1: Means (± standard deviations) of computational times taken by each method forthe respective slice images in Fig. 3.7.
Human data computational time (seconds)
L

R

LSC

0.161 ± 0.033

0.156 ± 0.016

MRTM2

1.045 ± 0.031

1.014 ± 0.031

LGA

0.073 ± 0.024

0.083 ± 0.009

The two columns indicated by L and R denote the left and right panel images, respectively.
patients’ images.

工学院大学

53

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

'

$

(a)

&

(b)

%

Figure 3.6: Comparison of the distributions of LSC errors (a) against OLS errors (b) for all
1024 simulated noisy tTACs for BPN D = 3.
工学院大学

54

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

'

$

&

%

Figure 3.7: Axial slice images of brain regions obtained by three methods, LSC, MRTM2, andconventional OLS-based LGA. These are BPN D images of Aβ quantification of 11 C-PiB brain
PET data of the two patients. The left panel shows the images of the Aβ-negative patient, andthe right panel shows the images of the Aβ-positive patient.

工学院大学

55

KOGAKUIN UNIVERSITYPK SHIGWEDHA

3.2

LSC

Discussion

The results in the above analysis shows that LSC estimates of the BPN D emerged to be lessbiased than those of MRTM2 and conventional OLS-based LGA. In addition, it was also seenthat LSC estimates have smaller variations in comparison with those for MRTM2, providing abetter trade-off between bias and variability.

The underestimation which is observed in the simulated data for conventional OLS-based LGA
estimates is also reflected in the real data, and accordingly, the reduction of the underestimation
seen in the simulated data for the LSC method is also reflected in the real data. Furthermore,
the underestimations observed for conventional OLS-based LGA method in both simulation andhuman data studies are in consistency with previous studies [33, 52], and this strengthens thereliability of the results reported here.

Both LSC and OLS are based on the assumption that the errors in the measurements are normally distributed. However, this assumption may not always be satisfied. A visual inspectionof the distributions of the LSC error values and the OLS residuals —in Figs. 3.5 and 3.6—
showed that the LSC error values have a better normal distribution compared with that of the
OLS residuals. This can be attributed to the fact that LSC takes into account the varying
magnitude of errors at each data point as demonstrated in [14].

LSC considers the errors in both variables without altering the original LGA equation, thus
maintaining its simplicity and inducing less variability. In addition, LSC does not only consider
the errors associated with term CT (t), but also the errors associated with each LGA variable
as a whole. The effort to include more details in LSC has increased the number of terms in Eq.
(3.1), compared with the fundamental OLS regression in Eq. (2.49), which only comprises theresiduals of the response variable. This could be a factor in the slight increase in the variabilityof LSC-based estimates compared to conventional OLS-based LGA estimates.

MRTM2 rearranges the conventional LGA equation into a multilinear regression form so thatthe noted noisy term, CT (t), only appears in the response variable [9, 10]. This only considers
the errors associated with CT (t). However, the rest of the terms in the LGA equation are notentirely noise-free; they are less noisy compared with CT (t). Not considering the noise in the
工学院大学

56

KOGAKUIN UNIVERSITYPK SHIGWEDHA

LSC

rest of the terms could lead to the bias in the MRTM2 estimates.

Based on the computational times of the three methods for obtaining slice images, the conventional LGA has the highest computational efficiency, followed by LSC, and then MRTM2.
The computation time for LSC would depend on the number of iterations, which in turn depends on the iteration condition, with the computation time increasing with the number of
iterations. All LSC images in Fig. 3.7 required at most three iterations. The iteration limit00
ing condition in this study, (αi0 − αi+1)/αi+1< 10−10 , was arbitrarily chosen. However, it was

found that by choosing a stricter condition, say 10−11 , the number of iterations did not increase.

Although the weighting formulas in Eq. (3.66) are inverse squares of the measurements, they
will not cause the usual problem of overproportional weighting of smaller measurements. Thisis because they are not the direct weighting coefficients. The final weighting coefficients are appearing in Eq. (3.51) as Wi Bi . Because they appear in both the numerator and denominator,
the overprortional weighting effect can be canceled out.

The proposed method, LSC-based LGA, significantly reduced the bias in the estimates of BPN D ,up to 12 and 28% difference compared with MRTM2 and conventional OLS-based LGA, respectively. Considering the conventional OLS-based LGA as the standard method, LSC reduced
the bias at a slight expense of increased variation. However, LSC caused smaller variations andrequired shorter computation times than MRTM2. Thus, LSC outperformed MRTM2 in all
aspects (bias, variability, and computational efficiency).

工学院大学

57

KOGAKUIN UNIVERSITY

Chapter 4

Tissue time-activity curves denoisingtechniques for improving image
contrast
4.1

Principal component analysis

PCA is a dimensionality reduction technique that is mostly used for feature extraction and noise
filtering. The fundamental idea of PCA is to represent the data in a reduced dimension while
preserving as much variance as possible, i.e., reducing a large data set of variables to a smaller
set that retains most of the important information in the original data set. PCA transforms
variables into new sets of variables that are linear functions of the original variables [15, 53, 54].
The new set of variables are uncorrelated and are defined by a set of orthogonal basis vectors
and PCs that optimally describe the variance in the data. That is, the first PC describes thelargest variance in the data. The second PC is orthogonal to the first PC and describes thelargest of the remaining variance. It goes on until the last PC – orthogonal to all earlier PCs –that describes the least variance than all other PCs. The orthogonal basis vectors and the corresponding variance are found by solving for the eigenvalues and eigenvectors of the covariancematrix of the data. The eigenvectors are the orthogonal basis vectors, and the eigenvalues are
the corresponding variances.

In dimension reduction, one keeps only the projections along the dimensions of the first few
largest PCs. This gives a reduced dimensional projection that maintains most of the important

58

PK SHIGWEDHA

PCA AND CORRCA

information of the original variables, allowing for a simpler interpretation.

In noise filtering application, one attempts to recover the original variables using the lowerdimensional projections. This results in denoised estimates of the original variables. The idea
is that the larger PCs correspond to the dimensions with most of the signals of interest and lessnoise. On the other hand, smaller PCs correspond to the dimensions with much of the unwanted
noise and less of the signal of interest. Thus, reconstructing the data using only the few largest
PCs, preserves the signals of interest whilst throwing away the unwanted noise signals.

Let us denote the PET brain volume data C to be made up of R tTACs with p frames. Then
we can write, c1 (t1 ) c1 (t2 ) . . .
 c2 (t1 ) c2 (t2 ) . . .

C= .
..
..
 ..
.
.


cR (t1 ) cR (t2 ) . . .
h= C(t1 ) C(t2 ) . . .


c1 (tp ) 
c2 (tp ) 

.. 
. 


cR (tp )C(tp )

(4.1)

R×p

i

Applying PCA to the data matrix C, it finds a set of p × p vectors, V ∈ Rp×p . These vectors
are simply the eigenvectors of the covariance matrix of the mean-centered data, C0 . Using only
the first k (k ≤ p) eigenvectors, a new set of denoised original variables can be estimated by,Â = (C0 W)W0 + C̄,(4.2)

where, W = V(:, 1 : k) ∈ Rp×k , is the set of the first k eigenvectors. C̄ ∈ R1×p is the mean vectorof C. The eigenvectors define the directions of the PCs, Z ∈ RR×p , which can be expressed as,
Z = (C0 W).

(4.3)

Therefore, using the first few k eigenvectors can also be referred to using the first few k PCs.

In this study, the number of PCs, k, that were retained to estimate the denoised data were suchthat specified amount of percentage variance could be realized from them. Specifically, the two
methods LSC–PCA and OLS–PCA were evaluated at two levels of percentage variance, 95%
and 97%.

工学院大学

59

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

Mathematics of PCA
Given our data matrix C ∈ RR×p , the aim is to find a linear combination of the columns
of matrix C with maximum variance [15]. Let us denote these linear combinations by Z,
i.e., the first PC Z1 is defined as,
Z1 = Ca1

=

p
X

(4.4)
ai1 C(ti ),i=1hwhere a1 is a vector of constants coefficients, a11 a21 . . .

ap1i∈ R1×p .

In general, the k th PC can be expressed as,
Zk = Cak

=

p
X

(4.5)
aik C(ti ).

i=1

Equation (4.4) is subject to ”var(Z1 ) is maximum”, and therefore Eq. (4.5) is subject
to ”var(Zk ) is maximum”. Further, the PCs are such that they are uncorrelated witheach other (i.e., cov(Zk , Zl ) = 0 | k > l ≥ 1), and ak has a unit length (i.e., aTk ak = 1)
[53, 15].

The variance of the first PC is given by,var(Z1 ) = hZ12 i − hZ1 i2=pX

ai1 aj1 hC(ti )C(tj )i −

i,j=1

p
X

ai1 aj1 hC(ti )ihC(tj )i

i,j=1

(4.6)
=

p
X

ai1 aj1 Sij

i,j=1= aT1 Sa1 ,where S is the sample covariance matrix, and Sij is the covariance of C(ti ) and C(tj ),respectively given by,
S=

工学院大学

1
CT C0
R−1 0

60

(4.7)

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

and
Sij = hC(ti )C(tj )i − hC(ti )ihC(tj )i = σC(ti )C(tj ) .

(4.8)

Identifying the linear combination with maximum variance, var(Z1 ), is equivalent tomaximizing the equation below [15],aT1 Sa1 − λ1 (aT1 a1 − 1),(4.9)

where λ1 is a Lagrange multiplier. Differentiating Eq. (4.9) with respect to a1 andequating to zero gives,
Sa1 − λ1 a1 = 0(4.10)(S − λ1 Ip )a1 = 0=⇒ Sa1 = λ1 a1 ,where Ip is a p × p identity matrix. Equation (4.10) implies that a1 is a unit-norm
eigenvector, whereas λ1 is the corresponding eigenvalue of the covariance matrix S [15].
λ1 is the largest eigenvalue of S, and thus the first pricipal component Z1 accounts forthe largest of variances in the data [53].

hThe second vector of coefficients, a2 = a12 a22 . . .

i1×p
ap2 ∈ R , is found by max-

imizing the variance, var(Z2 ), of the of the second PC, Z2 —under the constraints
cov(Z2 , Z1 ) = 0 and aT2 a2 = 1. The variance of Z2 , and the covariance of Z1 and Z2
are respectively given by,var(Z2 ) = aT2 Sa2 ,(4.11)

andcov(Z1 , Z2 ) = aT1 Sa2(4.12)= λ1 aT1 a2(using Eq. (4.10)).

Maximizing var(Z2 ) under the aforementioned constraints equates to [15],aT2 Sa2 − λ2 (aT2 a2 − 1) − φaT2 a1 ,
工学院大学

61

(4.13)KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

where λ2 and φ are Langrage multipliers. Differentiating Eq. (4.13) with respect to a2
gives,
Sa2 − λ2 a2 − φaT1 a1 = 0.

(4.14)

Multiplying through Eq. (4.14) on the left by aT1 we obtain,aT1 Sa2 − λ2 aT1 a2 − φaT1 a1 = 0.

(4.15)

From Eq. (4.15) we notice that aT1 Sa2 = λ2 aT1 a2 = 0 and aT1 a1 = 1, and thus φ = 0.
Equation (4.15) can then be written as,
Sa2 − λ2 a2 = 0(4.16)=⇒ (S − λ2 Ip )a2 = 0.
Equation (4.16) implies that a2 is also an eigenvector of the covariance matrix S,
whereas λ2 is the corresponding eigenvalue. λ2 is the second largest eigenvalue of S,
and thus the second pricipal component, Z2 , accounts for the second largest of variancesin the data [53].

Generally, the variance of the k th PC Zk is given by,var(Zk ) = aTk Sak ,(4.17)and maximizing var(Zk ) under the aforementioned constraints, ({cov(Zk , Zl ) = 0 | k >l ≥ 1}, and {aTk ak = 1 | k = 1, 2, · · · , p}), we get the k th eigenvector (ak ) and thecorresponding eigenvalue (λk ). In the same way as the eigenvectors and eigenvalues forthe first and second PC, λk is the k th largest eigenvalue of the covariance matrix S, andthe k th PC, var(Zk ), accounts for the k th largest fraction of the variance in the data [53].

hIn total, with our data matrix, C = C(t1 ) C(t2 ) . . .
hi
PCs as, Z = Z1 Z2 . . . Zp ∈ R × p, such that,
Z = CV,

工学院大学

62

i
C(tp ) ∈ R × p, we define the(4.18)KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

where,
a11 a21 . . .
a12 a22 . . .

V= .
..
..
 ..
.
.


a1p a2p . . .


ap1 

ap2 

.. 
. 


app

,

(4.19)

p×p

is an orthogonal matrix, whose columns are the eigenvectors of the of S.

The original matrix can be re-estimated from the PCs by,Ĉ = ZV0(4.20)= (CV)V0
If the PCs are calculated using the mean-centred data, C0 , (i.e., Z = C0 V), then Eq.
(4.20) becomes,Ĉ = (C0 V)V0 + C̄.

(4.21)

Employing only the first few eigenvectors, W = V(:, 1 : k), Eq. (4.20) then matches
Eq. (4.2) as below,
Ĉ = (C0 W)W0 + C̄
(4.22)= Â,and this gives the denoised lower dimensional estimates of the original data set C.

4.1.1

Simulation studies

Simulation data design
In addition to the CFN radiotracer simulations used in section 3.1.1, a simulation of 11 C-PiB
was also performed in this section. 11 C-PiB is a reversible Aβ-binding radiotracer, and it is
the radiotracer of the clinical data in both Chapter 3 and the current Chapter (Chapter 4).
Therefore, for simulation studies, the methods are evaluated for two radiotracers, CFN and
11

C-PiB, in this Chapter. In this Chapter, the methods are compared in terms of DVR, insteadof BPN D . It is noted earlier in Eq. (2.41) that BPN D is related to DVR by BPN D = DVR − 1.

The set of kinetic parameters — exchange rate constants K1 , k2–4 and the non-displaceable
工学院大学

63

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

distribution volume (VN D ) — below was used to make the first simulation data set (CFN
radiotracer). These kinetic parameters were adopted from a previous study [52]. These are also
the same range of parameters adopted by other studies such as [16].
• 11 values of DVR were set in the range of [0.0, 4.0], covering a range of [0.0, 0.35] for k3 .
• [K1 k4 ] = [0.1835 mL cm−1 min−1• [k2 k3 ] = [K1 /VN D• VN D

0.115 min−1 ].

(DVR − 1) · k4 ].

= 1.59 mL cm−1 .

The above parameters are similar to the ones in Chapter 3, except for that now there is DVR
instead of BPN D . Simulated data were made for a 60 min, 27-frame scan (6 × 0.1667 min,3 × 0.5 min, 5 × 1 min, 5 × 2.5 min, and 8 × 5 min). Two-tissue compartment PET data were
simulated using a clinically measured plasma time-activity curve. The reference region wasformed with a minimal noise level, and DVR = 1. As per the physiological assumptions of thereference region, the delivery (K1 ) and clearance (k2 ) rates were set identical to those of thetarget tissues.

For 11 C-PiB simulations, the kinetic parameters, K1 , k2–4 and VN D , below were used to simulate
two-tissue compartmental model PET data. These parameters covers a range of parameters in
one of the pioneering studies of 11 C-PiB [26].
• 11 values of DVR were set in the range of [0.0, 4.0], covering a range of [0.0, 0.045] fork3 .
• [K1 k2 ] = [0.262 mL cm−1 min−1• [k3 k4 ] = [(DVR − 1) · k4
• VN D

0.121 min−1 ].

0.015].

= K1 /k2 .

The scan period and frame for 11 C-PiB simulation was made similar to that of clinical 11 CPiB PET images used in this study. These were 70 min, 25-frame scan (6 × 0.1667 min,3 × 0.3333 min, 2 × 1 min, 2 × 3 min, and 12 × 5 min). The reference region was formed in asimilar manner as for CFN; minimal noise level; DVR = 1; and K1 and k2 identical to those of
the target tissues.

工学院大学

64

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

For each of the two radiotracers, a noise-free tTAC was formed for each DVR value. Statisticalnoise was added to the noise-free tTACs to form noisy tTACs. For each noise-free tTAC, 1024
noisy tTACs were simulated as a slice of 32 × 32 voxels. Four methods, LSC–PCA, OLS–PCA,LSC and conventional OLS-based LGA, were then used to estimate the DVR values from thenoisy tTACs, and the results were compared. The true value for the dissociation rate of thereference region (k2R ) used to construct the simulation data was used for all methods. The noisemodel used here is the same as descibed in section 3.1.1.
Results
Figure 4.1 presents the simulated tTACs, mimicking the two radiotracers, 4.1(a) CFN and4.1(b) 11 C-PiB. Noise-free tTACs are represented by the thick curves, whereas the thin curves
represent the noisy tTACs. For each noise-free curve, out of the 1024 simulated noisy curves,
only one is shown in each of the subfigures in Fig. 4.1. In both Figs. 4.1(a) and 4.1(b), thelowermost thick curve denotes the reference region, and the upper 11 curves denote the 11 DVRvalues. DVR estimations from 11 C-PiB appeared to be sensitive to noise in comparison to CFN,
and therefore, 11 C-PiB tTACs were simulated with lower noise levels than CFN, as it can beseen in Fig. 4.1. Lower radioactivity measurements in 11 C-PiB could be a contributing factor
to the raised sensitivity to noise.

Figure 4.2 shows examples of noisy tTACs, and their corresponding PCA-denoised tTACs.
Three denoised tTACs are shown for each radiotracer, and they represent PCA-denoising witha number of PCs that explains 80%, 95% and 97% of variance in the data, respectively. For
CFN, the tTACs denoised with 95% and 97% variances are still noisy to a noticeable extent.
As the number of PCs is lowered to explain only 80% variance, it can be seen that the tTACs
are further denoised. One still needs to be careful not to lower the number of PCs to a greater
extent. The tTACs would seem well-denoised, however spatial information will be lost and thiscannot be inferred from the tTACs. For 11 C-PiB, the tTACs are denoised to a greater extent
in comparison to CFN. This could be, in part, due to that they consists lower noise levels.

Figure 4.3 compares the percentage recovered DVR in the averages of the values estimatedfrom the simulated noisy tTACs using four methods, LSC–PCA, LSC, OLS–PCA and OLS. In
other words, the difference from 100% denotes the percentage bias. For PCA denoising, CFN

工学院大学

65

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

'

$
104

20

104

10
8

activity [Bq/mL]

activity [Bq/mL]

15

10

6
4
2

5
0
0

-2
0

10

20

&

30

40

50

60

0

10

20

30

40

time [min]time [min](a) CFN(b) 11 C-PiB

50

60

70

%

Figure 4.1: Simulated tTACs for the radiotracers, (a) CFN and (b) 11 C-PiB. In each figure,the thick and thin curves represent the noise-free and noisy tTACs, respectively. The upper 11
thick curves correspond to the 11 DVR values, whereas the 12th and lowest curve correspondto the time-activity curve of the reference region.

$

'
104

8

104

10

tTACS
Noise-free
Noisy
PCA-80%
PCA-95%
PCA-97%7
8

activity [Bq/mL]

activity [Bq/mL]

6
5
4
tTACS
Noise-free
Noisy
PCA-80%
PCA-95%
PCA-97%

3
2
1

6

4

2

0

0
0

&

10

20

30

40

50

60

0

10

20

30

40

time [min]time [min](a) CFN(b) 11 C-PiB

50

60

70

%

Figure 4.2: Simulated noisy tTACs and their corresponding PCA-denoised tTACs for the two
radiotacers, (a) CFN and (b) 11 C-PiB. Three denoised tTACs are shown for each radiotracer.
The three denoised tTACs correspond to PCA with a number of PCs that explains 80%, 95%
and 97% variance, respectively.

工学院大学

66

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

tTACs were denoised with 95% and 97% variances for the results in Figs. 4.3(a) and 4.3(b),respectively. The same is true for 11 C-PiB in Figs. 4.3(c) and 4.3(d). As expected, LSC–PCA
estimates have smaller error bars than LSC estimates. This is especially visible in the figuresfor 11 C-PiB, Figs. 4.3(c) and 4.3(d). This demonstrates PCA’s effectiveness of reducing thevariance. PCA’s variance reduction can also be seen for OLS–PCA estimates. With regards to
Figs. 4.3(c) and 4.3(d), it can be seen that OLS–PCA is greatly affected by the choice of thepercentage variance that should be accounted for in PCA denoising. With only 2% differenceof the explained variance, the difference in OLS–PCA estimates in Figs. 4.3(c) and 4.3(d) is
clearly visible.

Table 4.1 shows some of the average values of the estimated DVR for CFN raditracer, with
tTACs denoised with 95% variance PCs for PCA methods. The data in Tab. 4.1 correspondto some of the values plotted in Fig. 4.3(a). Table 4.1 also confirms reduced variance in both
PCA-based methods, and also the difference in bias reduction of these two PCA-based methods,
LSC–PCA and OLS–PCA, i.e., smaller bias for LSC–PCA compared to OLS–PCA.

Table 4.1: Averege DVR values (and standard deviations in brackets) estimated from the noisytTACs by the four methods, LSC–PCA, LSC, OLS–PCA and OLS, for CFN radiotracer. These
data correspond to Fig. 4.3(a), in which PCA methods are denoised with 95% variance PCs.
True DVR

1.250

2.350

2.900

3.450

4.00

Estimated DVR
LSC-PCA

1.249 (0.063)

2.356 (0.167)

2.895 (0.246)

3.456 (0.324)

3.992 (0.394)

LSC

1.249 (0.068)

2.344 (0.190)

2.881 (0.288)

3.452 (0.372)

4.022 (0.518)

OLS-PCA

1.212 (0.061)

2.233 (0.152)

2.719 (0.219)

3.212 (0.277)

3.656 (0.338)

OLS

1.211 (0.065)

2.212 (0.165)

2.683 (0.227)

3.169 (0.296)

3.616 (0.395)

工学院大学

67

KOGAKUIN UNIVERSITYPK SHIGWEDHA

4.1.2

PCA AND CORRCA

Clinical 11 C-PiB PET data studies

Study design
A cohort of 12 (11 Aβ-negative and 1 Aβ-positive) subjects was included in this study. The
subjects underwent 70 min 11 C-PiB PET scan of 25 frames (6 × 0.1667 min, 3 × 0.3333 min,2 × 1 min, 2 × 3 min, and 12 × 5 min). The scan procedure was carried out according to
the Alzheimer’s Disease Neuroimaging Initiative (ADNI) 11 C-PiB standards. The study wasapproved by the Ethics Committee of Kindai University Hospital, and written informed consent
was obtained from all participants. Obtained PET data were voxels of size 2.1 × 2.1 × 3.4 mm3 ,in arrays of 128 voxels × 128 voxels × 47 slices. DVR parametric images were then generated
from the PET data using the LSC–PCA, OLS–PCA, LSC, and OLS methods.

In addition to visual inspections, the parametric images were compared numerically in termsof the contrast between the four main gray matter cortices (frontal, temporal, occipital, andparietal) and white matter (corona radiata). The contrast index was calculated by the equation,C=µ W − µG,
σW + σG(4.23)

as in [55]. µG and µW denote the mean DVR of the gray and white matter regions, respectively,
whereas σG and σW denote the standard deviations of DVR in the gray and white matterregions, respectively. All calculations were performed in MATLAB (The MathWorks, Inc.,
Natick, MA, USA).
Results
Figure 4.4 shows a comparison of the medians of contrast between gray matter cortices andwhite matter for DVR images obtained by the four methods (LSC–PCA, OLS–PCA, LSC and
OLS), calculated for 11 Aβ-negative subjects. Statistical differences of these values are expressesd in Tab. 4.2. In the left figure, 4.4(a), the data were denoised with PCs that explained 95%
variance for the PCA-based methods. The four gray matter regions show a similar trend, withthe contrast of LSC–PCA and OLS–PCA images being the highest, as expected. As similar
results were observed for LSC–PCA and OLS–PCA in the simulation results for 11 C-PiB datadenoised with 95% of variance, here similar contrast is also observed between LSC–PCA and
OLS–PCA. The figure on the right, 4.4(b), shows contrast values of a single subject, and forwhich the tTACs were denoised with 97% of variance for the PCA methods. Numerical values
工学院大学

68

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

for this figure are given in Tab. 4.3. Similar to simulation results, here it is observed again that
OLS–PCA can perform poorly in comparison to LSC–PCA if more PCs are used. This serve to
inform that, more precaution is required for OLS–PCA to balance between bias and variancereduction. On the other hand, for LSC–PCA, one needs only to consider throwing away a few
PCs that could help with variance reduction, and not worry about the bias.

The Friedman Test Post-Hoc multiple comparison p-values of the difference between the contrasts of the images obtained using the four methods are shown in Tab. 4.2. The contrasts
being compared in Tab. 4.2 are those used to plot Fig. 4.4(a). According to these statisticalvalues, the contrasts of both LSC–PCA and OLS–PCA images are significantly different fromthose of LSC and OLS.

Table 4.2: Friedman test and Post-Hoc multiple comparisons of the p-values of the contrast.
Tabulated numerical values correspond to the p-values of the difference between the two methods
in the corresponding rows and columns.
Friedman Test (p value = 0.001)
Post-Hoc Multiple Comparisons p-values

LSC–PCAOLS-PCA

LSC

OLS

0.656

0.005

< 0.0010.025< 0.001

OLS–PCA
LSC

0.086

Figure 4.5 shows axial slices of the DVR parametric images obtained using the four methods.
PCA-based methods were denoised with the number of PCs that explained 95% variance forthese images. The two images in the upper row are from two Aβ-negative subjects, while thebottom row is from the Aβ-positive subject. In the figure, noise components in the images
estimated by the LSC method can be seen appearing as random noisy pixels/regions with high
DVR values throughout the slice (see the black arrows in Fig. 4.5). Images obtained by the
LSC–PCA method exhibit a rather fair distribution of the DVR values with no sharp increases

工学院大学

69

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

Table 4.3: Contrast values for a single Aβ-negative subject, which are plotted in Fig. 4.4(b).
Contrast
FrontalTemporal

Occipital

Parietal

LSC-PCA

1.162

1.030

1.110

1.049

OLS-PCA

1.030

0.939

0.962

0.847

LSC

1.007

0.915

0.919

0.789

OLS

0.905

0.831

0.813

0.624

in the pixels/regions. The images obtained by OLS–PCA appear similar to those for LSC–PCA.
This was expected since similar trends were observed between LSC–PCA and OLS–PCA for
DVR estimates in simulation results and also for contrast values, for the tTACs denoised with
95% variance. The LSC images thus appear slightly thick (in comparison to the LSC–PCA images), which can hinder the contrast between the brain structures. For the OLS-based images,
the lowest DVR estimates were observed as expected.

Figure 4.6 shows a single slice of DVR image for an Aβ-negative subject from which the contrasts
in Tab. 4.3 were calculated, for each of the four methods. The difference between LSC–PCA
and OLS–PCA can again be seen in this figure, with the OLS–PCA image being biased incomparison to that of LSC–PCA.

工学院大学

70

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

110

115
110

105
105
100

100
95

95

90
85
80

LSC-PCALSC
OLS-PCA
OLS

1.5

2

90

85
2.5

3

3.5

4

1.5

(a) CFN (PCA: 95% variance)
100

95

95

90

90

85

85

80

80

75

65

2

2.5

3

3.5

4

(b) CFN (PCA: 97% variance)

100

70

LSC-PCALSC
OLS-PCA
OLS

75
LSC-PCALSC
OLS-PCA
OLS

1.5

2

70
65
2.5

3

3.5

4

LSC-PCALSC
OLS-PCAOLS

1.5(c) 11 C-PiB (PCA: 95% variance)

2

2.5

3

3.5

4

(d) 11 C-PiB (PCA: 97% variance)

Figure 4.3: Percentage recovered DVR of the average (N = 1024) DVR values estimated fromthe noisy tTACs by the four methods, LSC–PCA, LSC, OLS–PCA and OLS. The error bars
denote the standard deviations. In the upper panel, (a) and (b), CFN tTACs were denoisedwith PCA at 95% and 97% variance respectively. The same applies for 11 C-PiB in the lower
panel, (c) and (d).

工学院大学

71

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA'$
Grey matter / white matter contrast

2

2

LSC-PCAOLS-PCA
LSC
OLS

1.8
1.6
1.4

LSC-PCALSC
OLS-PCA
OLS

1.8
1.6
1.4
1.2

Contrast

1.2

Contrast

Grey matter / white matter contrast

1
0.8

1
0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2
FrontalTemporal

Occipital

Parietal

Frontal(a) Clinical 11 C-PiB (11 Aβ-negative&

Temporal

Occipital

Parietal(b) Clinical 11 C-PiB (single Aβ-negative

subjects)

subject)%

Figure 4.4: (a) Box plots of the DVR contrast between gray matter regions, frontal, temporal,
occipital and parietal cortices, and white matter region, corona radiata, calculated for the 11
Aβ-negative subjects. For the PCA methods, the tTACs were denoised with 95% variance. The
circles in the box plots denote the means, whereas the horizontal lines denote the medians. (b)Plots of the DVR contrast of a single subject. For these plots, for the PCA methods, the brain
volume was denoised with 97% of variance.

工学院大学

72

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

2.5

2

1.5

1

0.5

0

Figure 4.5: DVR parametric images as obtained by the four methods, LSC–PCA, OLS–PCA,LSC and OLS. The two upper rows are for Aβ-negative subjects, whereas the lowermost row is
for the Aβ-positive subject. For PCA-based methods, the tTACs were denoised with PCs that
explained 95% of variance.

工学院大学

73

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA'$
LSC-PCAOLS-PCA

5

4

3
LSC

OLS

2

1

0

&

%

Figure 4.6: DVR parametric images as obtained by the four methods, LSC–PCA, OLS–PCA,LSC and OLS. In this case, the tTACs were denoised with PCs that explained 97% of variancefor the PCA-based methods.

工学院大学

74

KOGAKUIN UNIVERSITYPK SHIGWEDHA

4.2

PCA AND CORRCA

Correlated component analysis

Like PCA, CorrCA is also a feature extraction method. CorrCA operates by identifying components that are maximally correlated between repetitions in multivariate data [19]. Specifically,
CorrCA maximizes the ratio of between-repetition to within-repetition covariance. In the context of our application to PET data, our repetitions are along the slices dimension. CorrCA will
therefore maximize the ratio of between-slices to within-slices covariance. This ratio is generally
referred to as inter-subject correlation (ISC) [19]. Specifically in our context, it translates tointer-slice correlation.

Consider a set of dynamic PET brain volume data arranged as an array of size q × p × M ,where M denotes the number of slices, q denotes the number of voxels in a slice, and p denotes
the number of time points (frame numbers). CorrCA identifies directions in the p-dimensional
space along which the tTACs maximally correlate between M slices, with correlation measured
across q voxels.

Let us put it into perspective in comparison to PCA applied to the same data. For the PCA
case, the brain volume will be represented as an R × p array, where R is a product of q and
M , and denotes the total number of voxels in the whole brain volume. PCA returns a
set of p-dimensional vectors which successively capture the variance in the data in a descendingorder. Similarly, CorrCA returns a set p-dimensional vectors which successively capture the ISC
in a descending order. CorrCA also formulates into an eigenvalues and eigenvectors problem,from which the p-dimensional orthogonal vectors (and correlated components) are found as the
eigenvectors, and ISCs as the eigenvalues.

Maximizing the between-slices to within-slices covariances maximizes the mean-over-varianceacross slices, which has been asserted to define a signal-to-noise ratio [19]. CorrCA capitalizeon the ability to simultaneously operate ”within individual slices (within-slices covariances)”and ”across all slices (between-slices covariances)”, without arranging the slices in one plane.
That means operating directly on the q × p × M array, unlike for PCA where the slices arearranged in one plane side-by-side to form the R × p array. This means that CorrCA is a dual
operation, exploiting through all slices whilst also treating each slice individually. This can beseen as advantageous over PCA, which means CorrCA can possibly provide a better noise filter
工学院大学

75

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

method over PCA. This therefore forms the motivation to introduce CorrCA to PET parametric
imaging and assess the performance of CorrCA in comparison to PCA. For CorrCA only twocorrelated components were retained to denoise the tTACs.

For a brief mathematical outline of CorrCA as presented in [19], let the observed noisy tTACs
be xji = [cji (t1 ), cji (t2 ), · · · , cji (tp )]T , wherei = 1, · · · , q and j = 1, · · · , M . (T denotes atranspose). Thus, xji ∈ Rp×1 is the ith voxel in the jth slice. This slice can be represented as,
xji (t) = Ayij + εji ,(4.24)

where A ∈ Rp×k (k ≤ p) is a projection matrix to be constructed by CorrCA; yij ∈ Rp×1 is thecoefficient vector; and εji ∈ Rp×1 is the residuals vector.

The jth slice can be written as,jc1 (t1 )cj1 (t2 )...
 jc (t1 ) cj (t2 ) . . .
2
 2
Cj =  .
..
..
 ..
.
.


cjp (t1 ) cjp (t2 ) . . .



cj1 (tp )


cj2 (tp )

.. 
. 


cjp (tp ),(4.25)

q×p

and the whole brain data as,h
C = C1 C2 . . .

CM

i

.

(4.26)

q×p×M

The objective of CorrCA is to find a linear combination of the p measurements (yij ), defined as,yij = vT xji ,(4.27)

such that the correlation across M slices is maximized —where v ∈ Rp×1 is a projection vector.
The ISC (ρ) is defined as the ratio of the between-subject covariance, rB , and the within-subject
covariance, rW , which can be expressed as,
ρ=

1
rB
.
M − 1 rW

(4.28)The term 1/(M − 1) is such that ρ ≤ 1. In our context, the between-slices covariance can beinterpreted as a sum over all pairs of slices, and within-slices covariance can be interpreted as
工学院大学

76

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

a sum over all slices, as in the equations below,
rB =

M
X

M
X

q
X
(yij − ȳj )(yik − ȳk )j=1 k=1, k6=j i=1(4.29)
rW =

q
M X
X

(yij − ȳj )(yij − ȳj ),j=1 i=1

where ȳj = 1q

Pqji=1 yi is the sample-mean of the projected data for the jth slice. Substituting

Eq. (4.27) into Eq. (4.29) gives,
rB = vT RB(4.30)
rW = vT RW ,where RB and RW are the between-slice and within-slice covariance matrices of xji , respectively
given by,
RB =

q X
M
X

M
X

(xji − x̄j∗ )(xki − x̄k∗ )T(4.31)

i=1 j=1 k=1, k6=j

RW =

q X
M
X

(xji − x̄j∗ )(xji − x̄j∗ )T .

(4.32)

i=1 j=1

Substituting Eq. (4.30) into Eq. (4.28) gives,ρ=vT RB v
1.
M − 1 vT RW v

(4.33)

To maximize ρ, we rearrange and equate Eq. (4.33) to zero and differentiate with respect tovT , i.e.,

∂
T
T
ρ
v
R
v
−
v
R
v
=
0
W
B
vT

(4.34)=⇒ ρ RW v = RB v.
Equation (4.34) means that v is an eigenvector of RW −1 RB , and ρ is the corresponding eigenvalue —under the assumption that RW is an invertible matrix. It is therefore that ISC is
maximized by projecting the data onto the eigenvector of RW −1 RB with the largest eigenvalue[19]. In a similar manner as for PCA, these eigenvector and eigenvalue can be labeled as v1
and ρ1 , denoting the first eigenvector and eigenvalue, respectively. Further in a similar way to
PCA, successive eigenvectors (v2 , · · · , vp ) and corresponding eigenvalues (ρ2 , · · · , ρp ) captures
工学院大学

77

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

the ISC in a decreasing order. To find all the eigenvectors and eigenvalues, maximal ISC issummed over all components. The components are required to be uncorrelated, which means
each component captures “different and independent” characteristics of the data [19].

In entirety, Eq. (4.33) can therefore be expressed as,
J(V) =

p
X

ρd =

d

p
X
vT RB vddd

vdT RW vd.

(4.35)

Thus, CorrCA maximize Eq. (4.35) with respect to V = [v1 , v2 , · · · , vp ], subject to vdT RW vc = 0
for c 6= d. The solution to Eq. (4.35), corresponding to Eq. (4.34) is,
RB V = RW VΛ.

(4.36)

Here, Λ ∈ Rp×p denotes a diagonal matrix consisting of the ISC values (ρd ) in the diagonals.
These ISC values are eigenvalues of the eigenvectors vd which defines the projections yd = vdT x.
The first projection, along the first eigenvector, provides the highest ISC across M slices. Successive projections are uncorrelated with each other and captures ISC in decreasing order, i.e.,the second projection correspond to the direction with second largest ISC and so forth, until
the last projection which corresponds to the direction with the smallest ISC. [19].

In matrix format, Eq. (4.24) can be written as,
X|p×q×M = (A|p×k ) (Y|k×q×M ) + ε|p×q×M .

(4.37)The least-squares estimate of A is,−1

Â = RW v(vT RW v).

(4.38)

Analogous to using only a few PCs in PCA, here, using only the first few k (k ≤ p) correlated
components with the highest ISC values in Eqs. (4.37) and (4.38), denoised tTACs can be
estimated as per Eq. (4.24) as,Ĉ = X̂ = ÂŶ.

(4.39)For mean-centred input data, C0 , we get that,Ĉ = ÂŶ + C̄,(4.40)

where C̄ is the mean of C.

工学院大学

78

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCAIn short, it can be inferred from Eqs. (4.34) and (4.36) that CorrCA performs PCA on the matrix of the ratio of between-slices to within-slices, RB R−1
W , whose the eigenvalues define the ISC.

CorrCA was computed using the MATLAB codes provided by Parra and colleagues [19]. The
code can be found at http://parralab.org/corrca.

4.2.1

Simulation studies

Simulation data design
The simulation parameters and conditions used for CorrCA assessment in this section are thesame with those used in section 4.1.1 for PCA studies. However, here a further step was taken toevaluate the methods at different noise levels. The magnitudes of the noise level were calculated
as the percentage of the ratios of the “standard deviation of the noisy tTACs” to the “mean of
the noise-free tTAC” for the portion of the time range used for the DVR estimation. For CorrCA
assessment, the focus is mostly on LSC-CorrCA against LSC-PCA. However, a brief analysis
of LSC–CorrCA and OLS–CorrCA, demonstrates (for OLS–CorrCA) the similar issue observedwith OLS–PCA that one needs to be extra careful when choosing the number components toretain to denoise the tTACs. The results reported for LSC-CorrCA are mainly for the datadenoised with only two components (correlated components). This is to assess if LSC-CorrCAwith only two correlated components can achieve better or similar results with LSC–PCA withmany PCs.
Results
Figure 4.7 shows a simulated noise-free tTAC, a corresponding noisy tTAC, and three denoisedtTACs, for the two radiotracers, CFN 4.7(a) and 11 C-PiB 4.7(b). Two of the denoised tTACs
were denoised by CorrCA with 2 and 5 components, as denoted by the subscripts in the legend.
For PCA, examples of tTACs denoised with many components are already shown in section
4.1.1 (Fig. 4.2). In this figure the PCA-denoised tTACs were denoised with 2 PCs. The tTACs
denoised by CorrCA with two components appears to to give the best estimate of the true
(noise-free) tTAC. The tTAC denoised by PCA with two PCs appear to overestimate the true
curve in Fig. 4.7(a) for CFN, but appears to underestimate the one for 11 C-PiB in 4.7(b).
Overall all three denoised tTACs appears to be considerably denoised. However, especially for
PCA, it is known that denoising the tTACs with only two PCs leads to information loss which
工学院大学

79

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

leads to low quality images. That means looking at the tTACs alone is not enough to conclude
whether a denoising method will lead to improved images.
'

$
104

10

10

8

tTACS
Noise-free
Noisy
CorrCA 2

8

activity [Bq/mL]

activity [Bq/mL]

104

6
tTACS
Noise-free
Noisy
CorrCA 2

4

2

CorrCA 5

6

PCA2

4

2

CorrCA 5
PCA2

0

0
0

&

10

20

30

40

50

60

0

10

20

30

40

time [min]time [min](a) CFN(b) 11 C-PiB

50

60

70

%

Figure 4.7: Simulated noisy tTACs and their corresponding CorrCA- and PCA-denoised tTACs
for the two radiotacers, (a) CFN and (b) 11 C-PiB. Three denoised tTACs are shown for eachradiotracer – two denoised by CorrCA and one denoised by PCA. The subscripts denote thenumber of components used to denoise the corresponding curves.

Figure 4.8 shows the percentage bias in the averages of DVR values estimated from the noisytTACs by four methods, LSC–CorrCA, LSC–PCA, LSC, and OLS. The percentage bias is thedifference from the 100% mark, i.e., the closer the plot is to 100% the less biased it is. Herethe focus is on LSC–CorrCA and LSC–PCA. It can be seen in both Figs. 4.8(a) and 4.8(b)that the variance is drastically reduced in the LSC–CorrCA estimates. This is mostly because
only 2 components were used to denoise the data for LSC–CorrCA, in comparison to many PCs(which could go up to 10) used to retain 95% variance for LSC–PCA. Estimates from tTACs
denoised by PCA with 2 PCs would also have smaller variances. It then remains to evaluate
the clinical images to assess which method can hold much information in a smaller number ofcomponents.

In Fig. 4.9, it is shown for both radiotracers that, when all the correlated components areused, CorrCA returns the original data, and therefore the DVR estimates of LSC–CorrCA are
identical to those of LSC. This is also a known characteristic of PCA [16] , i.e., when all PCs
工学院大学

80

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

'

$

115
100

110

95

105

90

100

85

95

80

90
LSC-CorrCA
LSC-PCA
LSC
OLS

85
80

1.5

2

&

75
70
2.5

3

3.5

4

LSC-CorrCA
LSC-PCA
LSC
OLS

1.5

2

2.5

3

(b) 11 C-PiB(a) CFN

3.5

4%

Figure 4.8: Percentage recovered DVR of the average (N = 1024) DVR values estimated fromthe noisy tTACs by four methods, LSC–CorrCA, LSC–PCA, LSC and OLS. The error bars
denote the standard deviations. In both figures, (a) and (b), the tTACs were denoised withtwo components for LSC–CorrCA. For LSC–PCA, the tTACs were denoised with the numberof PCs accounting for 95% of variance in the data.

are used PCA returns the original data.

Figure 4.10 shows the slices estimated from the noisy tTACs by four methods, LSC–CorrCA,
LSC–PCA, LSC and OLS. In this figure, it can be seen that the LSC–CorrCA estimates show
better estimates than other methoods, in agreements with the smaller variances observed inFig. 4.8. The noise-free slices are shown in Fig. 4.11.

Figure 4.12 compares the DVR values estimated by LSC–CorrCA and OLS–CorrCA, in which
the tTACs were denoised with 5 components. It is revealed here that if many components
are used, it raises the bias in the OLS–CorrCA estimates, similar to what was observed for
OLS–PCA. This further puts LSC to be the regression method of choice for LGA-based PETparametric imaging.

In Fig. 4.13, the four methods, LSC–CorrCA, LSC–PCA, LSC and OLS are compared at different noise levels. In Fig. 4.13(a), the estimates are plotted individually for each noise level.
In Fig. 4.13(b), the estimates at all noise levels are put together in a box plot. It is seen again
工学院大学

81

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

'

$

115

LSC-CorrCA (2 CCs)
LSC-CorrCA (All CCs)
LSC
OLS

110

100
95

105
90
100
85
95
80
90
75
85
70
80
1.5

&

2

2.5

3

3.5

4

LSC-CorrCA (2 CCs)
LSC-CorrCA (All CCs)
LSC
OLS

1.5

2

2.5

3

3.5

4

(b) 11 C-PiB(a) CFN%

Figure 4.9: Percentage recovered DVR of the average (N = 1024) DVR values estimatedfrom the noisy tTACs by four methods, LSC–CorrCA (with two correlated components), LSC–
CorrCA (with all correlated components), LSC and OLS. The error bars denote the standard
deviations. Figures (a) and (b) shows the estimates for CFN and 11 C-PiB radiotracers, respectively.

here how LSC–CorrCA has reduced variances.

In Fig. 4.13, however, the minimum number of PCs was set to three. This means that, if thenumber of PCs that accounts for 95% of variance are less than 3, 3 PCs are used, and if thenumber of PCs that accounts for 95% of variance are more than three, then that number of
PCs are used. Figure 4.14 is the equivalent of Fig. 4.13 without the restriction of minimum
number of PCs of 3. The effect of this condition can be seen for the first three lower noiselevels in Fig. 4.14(a), in which LSC–PCA shows an underestimation in comparison to all othermethods —including OLS-based LGA. The underestimation in LSC–PCA-based estimates isalso reflected in Fig. 4.14(b). For the first three lower noise levels, the number of PCs that
accounted for 95% of variance is one. This observation provides further insight into the notionthat using a few number of PCs leads to information loss. A similar trend was observed in[16], in which 1 PC-based estimates appeared slightly biased compared to those of 2 PCs. It istherefore important to always consider including more than 1 PC when denoising by PCA.

工学院大学

82

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA'$LSC-CorrCA

LSC-PCA

4
3.5
3
2.5

LSC

2

OLS

1.5
1
0.5
0

&

%

Figure 4.10: A display of the slices obtained by four methods, LSC-CorrCA, LSC-PCA, LSC and
OLS from simulated noisy tTACs for CFN radiotracer (DVR = 4). Each color block represents
a slice of 32 × 32 voxels, which makes a total of 1024 tTACs. The tTACs were denoised with 2
components for LSC–CorrCA, and with 95% variance for LSC–PCA.

4.2.2

Clinical 11 C-PiB PET data studies

Study design
The clinical data described and used for analysis in section 4.1.2 are the same data used for thissection. Under this section, the focus is completely on LSC–CorrCA vs LSC–PCA.
Results
Figure 4.15 compares the box plots of the contrast between gray and white matter for images
obtained by two methods, LSC–CorrCA and LSC–PCA. The tTACs were denoised with 2 components for LSC–CorrCA, and with 95% variance for LSC–PCA. The contrast was calculated
for 11 Aβ-negative images, as it was done in section 4.1.2 for Fig. 4.4(a). Some difference could
be seen in the mean (circle marks) and medians (horizontal lines) in the box plots of the two
工学院大学

83

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

'

$
1
0.8
0.6
0.4
0.2
0

&

%

Figure 4.11: A display of the simulated noise-free slices.

methods in Fig. 4.15. However the p value of their statistical differences across the four regions
as calculated by the Friedman Post-Hoc multiple comparison is 0.076, slightly over the desired
0.05 statistical significance mark. Numerical values for Fig. 4.15, mean values ans standard
deviations are shown in Tab. C.1.

Table 4.5 shows the computation time of four methods, LSC–CorrCA, LSC–PCA, LSC and
OLS. These data show another advantage aspect of the LSC–CorrCA method compared to
LSC–PCA, in which LSC–CorrCA takes shorter computation time compared to LSC–PCA. As
expected, both LSC–CorrCA and LSC–PCA took longer computation time compared to LSC.
This is because of the additional task of data denoising, by CorrCA and PCA, in addition toregression parameter estimation by LSC. The OLS approach required the shortest computation
time, as expected as well.

Figure 4.16 compares the images obtained by two methods, LSC–CorrCA and LSC–PCA, for two
subjects, Aβ-negative and Aβ-positive. For the Aβ-negative subject, the LSC–CorrCA image
appears to show higher DVR estimates. For the Aβ-positive subject, DVR estimates across the
slice appear about the same for LSC–CorrCA and LSC–PCA, with LSC–PCA image showing
a rather smoother distribution. For LSC–CorrCA images, it is worth noting that they were all

工学院大学

84

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

'

$

100
95
90
85
80
75
70

LSC-CorrCA
OLS-CorrCA
LSC
OLS

1.5

2

2.5

3

&

3.5

4
%

Figure 4.12: Comparison of the percentage recovered in the average (N = 1024) DVR values
estimated from the noisy tTACs by four methods, LSC–CorrCA, OLS–CorrCA, LSC and OLS.
For these results, the tTACs were denoised with 5 components for both LSC–CorrCA and
OLS–CorrCA.

estimated with 2 components, and an optimization for the number of components could makefor an improvement of the images. Overall, these results show that LSC–CorrCA with only 2components achieved about the same or better results with LSC–PCA with many components.
This means that CorrCA contains much information in a smaller number of components, incomparison to PCA. For PCA, the number of PCs that accounted for 95% variance in the 12
clinical data sets (11 Aβ-negative subjects and 1 Aβ-positive subject) used in this study were
all greater than 2, ranging from 4 to 8.

工学院大学

85

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA
$

'
4.6

7

4.4
6
4.2
5

4
3.8

4
3.6
3

3.4
3.2
0.0234 1.4256 2.8324 4.2173 5.6658

2

7.2995LSC-CorrCA

LSC-PCA

LSC

OLS

Noise Level

(a)

&

(b)

%

Figure 4.13: DVR estimates of CFN radiotracer. (a) A comparison of DVR (DVR = 4) estimation at different noise levels by four methods, LSC–CorrCA, LSC–PCA, LSC and OLS. (b) Box
plots of the estimated DVR (DVR = 4) values. These plots includes the values estimated at allnoise levels. The tTACs were denoised with two components for LSC–CorrCA. For LSC–PCA,
the tTACs were denoised with the number of PCs accounting for 95% of variance in the data(minimum number of PCs set to three).

'

$

4.6

6.5

4.4

6
5.5

4.2

5
4
4.5
3.8
4
3.6
3.5
3.4

3

3.2
0.0235 1.4082 2.8028 4.2165 5.6625

2.5

7.3074LSC-CorrCA

LSC-PCA

LSC

OLS

Noise Level

&

(a)

(b)

%

Figure 4.14: DVR estimates of CFN radoitracer. (a) A comparison of DVR (DVR = 4) estimation at different noise levels by four methods, LSC–CorrCA, LSC–PCA, LSC and OLS. (b) Box
plots of the estimated DVR (DVR = 4) values. These plots includes the values estimated at allnoise levels. The tTACs were denoised with two components for LSC–CorrCA. For LSC–PCA,
the tTACs were denoised with the number of PCs accounting for 95% of variance in the data.

工学院大学

86

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA'$

2LSC-CorrCA
LSC-PCA

1.8
1.6
1.4

Contrast

1.2
1
0.8
0.6
0.4
0.2
0
-0.2
FrontalTemporal

Occipital

Parietal&%

Figure 4.15: Box plots of the DVR contrast between gray matter regions, frontal, temporal,
occipital and parietal cortices, and white matter region, corona radiata, calculated for the 11
Aβ-negative subjects for the images obtained by two methods, LSC–CorrCA and LSC–PCA. For
LSC–CorrCA, the tTACs were denoised using two components. For the LSC–PCA, the tTACs
were denoised with 95% variance. The circles in the box plots denote the means, whereas thehorizontal lines denote the medians.

Table 4.4: Numerical values for the contrast between gray and white matter regions of theimages obtained by two methods, LSC–CorrCA and LSC–PCA.
Contrast
FrontalTemporal

Occipital

Parietal

LSC-CorrCA

1.236 ± 0.341

1.224 ± 0.360

1.136 ± 0.349

0.870 ± 0.266

LSC-PCA

1.105 ± 0.329

1.157 ± 0.333

1.056 ± 0.351

0.825 ± 0.296

工学院大学

87

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA'$LSC-CorrCA (A -negative)

LSC-PCA (A -negative)

2.5

2

1.5
LSC-CorrCA (A -positive)

LSC-PCA (A -positive)

1

0.5

0

&

%

Figure 4.16: DVR parametric images of a Aβ-negative and Aβ-positive subject, obtained by
two methods, LSC–CorrCA and LSC–PCA. The tTACs were denoised with two components
for LSC–CorrCA, and with 95% variance for LSC–PCA.

Table 4.5: Means (± standard deviations) of computation times taken by each of the four
methods, LSC–CorrCA, LSC–PCA, LSC, and OLS, to obtain a whole DVR image of size 128 ×
128 × 47 —calculated for the 11 Aβ-negative subjects and 1 Aβ-positive subject (i.e., N = 12).
Computation time (s)LSC–CorrCA

LSC–PCA

LSC

OLS

9.249 ± 1.811

10.678 ± 1.225

7.514 ± 0.879

2.570 ± 0.165

工学院大学

88

KOGAKUIN UNIVERSITYPK SHIGWEDHA

4.3

PCA AND CORRCA

Discussion

The results obtained in section 4.1 demonstrate that the LSC–PCA method improves the contrast of parametric images, in comparison with LSC. With LSC and PCA combination, LSC
reduces the bias, and PCA could help fend off the effects of the variations in the estimates,
resulting in parametric images with both minimal bias and variance. OLS–PCA was also foundto have reduced bias and improved contrast, but this was found to be greatly affected by the
number of PCs used to denoise the tTACs. For the CFN and 11 C-PiB simulations, it was seenthat the OLS–PCA estimates could get biased if one is not careful in determining the numberof PCs. In this case, one should then rely on LSC for bias reduction. This should therefore
highlight the importance of adopting a less-biased regression method, LSC, even when PCA
tTACs denoising method is applied. This is important especially when one is unsure of thenumber of PCs to use, which is usually the case in practical settings.

PCA has been shown to reduce bias using only a few PCs [16]. Using a minimal number of PCs
might result in information loss, leading to lower contrasts. This means that information is lostwith the number of PCs discarded, i.e. the fewer the number of PCs used the more information
is lost. The question of how many PCs to retain has no exact answer. It mostly depends on
the purpose of the study in question. For PET imaging, to maintain the spatial resolution of
the images, it is necessary to retain sufficient variance to reconstruct denoised tTACs. Thiscan always be set by determining the percentage variance to be accounted for in the denoised
data, and all the PCs that add up to that value can then be retained. For OLS–PCA it istherefore that using a few PCs leads to information being lost. On the other hand, using many
PCs reintroduces the noise in the data, and the estimates get biased. This is what LSC–PCA
is addressing; to allow many PCs to be used without risking bias.

The CorrCA denoising approach was applied to PET parametric imaging for the first time, andassessed against PCA. CorrCA is one of the many various approaches to canonical correlation
analysis (CCA) of multidata set, developed to identify linear combinations of signals that are
reliably reproduced across the data. In our context of a single PET image, CorrCA operates
by maximizing the ratio of the correlations measured within and across the slices in a PET
image. This is in contrast to PCA which focuses on the variance. Here, CorrCA is employed
as a temporal filter of tTACs, similar to the manner in which PCA is used. By using only a
工学院大学

89

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

few dimensions with the largest ISCs, one could get rid of the noise in the data along with theremaining components with smaller ISCs, whilst preserving the essential characteristics of thedata.

The simulation results showed that both LSC–CorrCA and LSC–PCA reduced the variance inthe DVR estimates, in comparison to LSC. This was indicated by smaller error bars in theplotted figures. LSC–CorrCA estimates showed the smallest error bars amongst all methods,
suggesting that it will give clinical images with improved contrasts. It was also shown thatOLS–CorrCA estimates get biased if more components are used, similar to OLS–PCA; thisinforms that LSC should be the preferred regression method to estimate the LGA slope.

The contrast of LSC–CorrCA images were observed to be higher than those for LSC–PCA images, as expected based on the simulation results. However, this difference was not statisticallysignificant in terms of the calculated p-value. These findings based on the calculated contrasts
are in agreement with the visuals of the displayed images, as not much of a difference could beseen between the images obtained by the two methods, LSC–CorrCA and LSC–PCA. It should
however be noted that LSC–CorrCA with only two components was able to achieve seemingly
better or similar results with LSC–PCA with many PCs. This could be that CorrCA capturemuch information in a less number of components in comparison to PCA. This will require further analysis to be established as such. As for computation time, LSC–CorrCA outperformed
LSC–PCA.

Another aspect of LSC–CorrCA is depicted in Fig. 4.17. Figure 4.17 shows scatters of DVR
estimates at different noise levels for four methods, LSC–CorrCA, LSC–PCA, LSC, and OLS.
In comparison to other methods, the variance of LSC–CorrCA estimates seems to be unaffectedby noise level, as the dispersion in the scattered DVR estimates does not show an increase withnoise level. This signifies a positive feature for LSC–CorrCA. Other methods show a clear increase in the dispersion of DVR estimates with respect to noise level. For LSC–CorrCA method
however, the scattered DVR estimates show some irregularities with respect to noise level, withsome estimates at lower noise levels exhibiting larger dispersion than those at higher noise levels.
This implies some instability in the LSC–CorrCA method. For CFN radiotracer, the effect of
this observation appears to be not of concern in comparison to LSC–PCA. For 11 C-PiB radio-

工学院大学

90

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

tracer, however, it appears these estimates with unexpected larger variances attains varianceslarger than those for LSC–PCA. The root cause of this instability resulting in random increases
in the dispersion of LSC–CorrCA estimates is unclear at this point, and it will be studied infuture researches.

In the approach proposed in this study, the number of correlated components to keep was determined with Â and Ŷ, which are used in Eqs. (4.39) and (4.40). As explained earlier, thisapproach is equivalent to performing PCA on the matrix RB R−1
W . In other words, this approach is as well equivalent to using a few, k, principal vectors in vT in Eqs. (4.27) and (4.38).
Therefore, the ISC values correspond to the variances explained by the eigenvectors of RB R−1
W
in this PCA operation.

An alternative method is to perform the truncation in RB and RW in Eq. (4.36) [19]. In
doing this, the resulting Â and Ŷ will also have reduced dimensions. This approach seemed tohave reduced the irregularities observed in Fig. 4.17. However, this approach required a larger
number of dimensions (correlated components) for better contrasts to be obtained in clinicalbrain images. Based on the analysis of real data contrast values, this approach also seemed tofit better with a method based on shuffle statistics [19] for determining the number of relevantcorrelated components. Shuffle statistics is used to test for nonlinearity of time series data. In
short, the idea of the shuffle statistics method is to generate surrogate data that are consistent
with the null hypothesis, but also similar to the original data in terms of time and spatial
correlation [19]. For CorrCA, in the context of our PET brain images, the null hypothesis
states that there is no shared one-dimensional subspace in which the data of at least two slices
are correlated. A shuffle statistics technique called random circular shift [19] is used to test
a significant ρ based on the consideration that there is no correlation between slices after arandom circular shift for each slice [19]. However, the number of correlated components thatgave better contrast values were smaller than those estimated by the shuffle statistics-based
approach —hence, overestimation. Overall, the approach used in this study (above paragraph)appears better than the approach described in this paragraph. However, a thorough analysis ofthese two approaches will be necessary in the future.

工学院大学

91

KOGAKUIN UNIVERSITYPK SHIGWEDHA

PCA AND CORRCA

'

$

7
8
6

7
6

5
5
4
4
3

3

2
2
LSC-CorrCA LSC-PCA

&

LSC

OLS

LSC-CorrCA LSC-PCA

LSC(b) 11 C-PiB(a) CFN

OLS%

Figure 4.17: Scatter of DVR estimates at 6 different noise levels (DVR = 4). (a) Scatter of
DVR values estimated by four methods, LSC–CorrCA, LSC–PCA, LSC and OLS, from noisy
CFN tTACs. (b) Scatter of DVR values estimated by four methods, LSC–CorrCA, LSC–PCA,LSC and OLS, from noisy 11 C-PiB tTACs. The tTACs were denoised with two components for
LSC–CorrCA. For LSC–PCA, the tTACs were denoised with the number of PCs accounting for
95% of variance in the data.

工学院大学

92

KOGAKUIN UNIVERSITY

Chapter 5

Summary and future directions
This study aimed at improving the LGA-based parametric estimates by reducing both thebias and variance. In the first part of the study (Chapter 3), an alternative linear regression
method, LSC, was employed, and its viability was assessed for BPN D estimation. The results
showed that LSC-based estimates are the least biased in comparison to those of MRTM2 andconventional OLS-based LGA. Specifically, LSC reduced the bias in the estimates of BPN D ,up to 12% and 28% difference compared with MRTM2 and conventional OLS-based LGA, respectively. In comparison to the standard OLS-based LGA, LSC reduced the bias at a slight
expense of increased variation. However, in total, LSC providied a better trade-off between bias
and variability, making it a promising tool for LGA-based PET parametric imaging.

Having noted the slightly increased variances in the LSC-based estimates in comparison with
OLS-based estimates, the second part of the study (Chapter 4) employed two tTACs denoisingtechniques. These are PCA and CorrCA. Since LSC-based estimates are already least-biased,
the two tTACs denoising methods were employed to reduce the variance in parametric estimates.
Both PCA and CorrCA operates by finding the most important dimensions. These approaches
were assessed in terms of DVR estimation. The results obtained demonstrate that both LSC–
CorrCA and LSC–PCA methods improved the contrast of the parametric images. This was
supported by a reduction in the variation of the estimates observed in the simulation results.
The contrast of LSC–CorrCA images were observed to be slightly higher than those for LSC–
PCA images, as expected based on the simulation results. On the other hand, OLS–CorrCA
and OLS–PCA were also found to improve the contrast of parametric images. However, morecaution is needed for the OLS-based method because if more components are used, the estimates93

PK SHIGWEDHA

SUMMARY

could get biased. This therefore demonstrates that it is important to always consider using LSC
as a regression method even when the data are denoised with CorrCA or PCA. Thus, this forms
the basis to encourage the adoption of LSC-based approaches, LSC–CorrCA and LSC–PCA, inwhich the complementary capabilities of LSC and PCA or CorrCA assures that both bias andvariance will be reduced, unlike in OLS–CorrCA or OLS–PCA in which one always risks bias.
These findings should be useful in PET parametric imaging for appropriate interpretation ofradiotracer binding in the brain. These are complementary approaches in which LSC mainly
serve to reduce the bias, and PCA and CorrCA are employed to minimise the variance in theestimates.

For LSC, a deeper look into weighting techniques could provide for improved estimates. The appropriateness of the weight function used in this study is demonstrated in section 3.1. However,they are estimates, and there is still a room for improvements. Improved weighting could allow
LSC to be a stand-alone approach, which means the data will no longer need to be pre-treatedwith methods like PCA or CorrCA prior to parameter estimation by LSC.

Simulation studies showed that the number of components required to adequately denoise thetTACs could vary depending on factors such as the radiotracer or noise level. Therefore, for
PCA and CorrCA, further studies could be to explore the methods to determine the numberof components to be retained to denoise the data. This is especially worth a while for CorrCAsince it is showing the smallest variations (in simulation data) and the highest contrasts (inclinical data); the aim will be to clearly reflect these observations in the displayed parametricimages. Investigation and comparison to other tTACs denoising methods would also make good
grounds for further studies.

工学院大学

94

KOGAKUIN UNIVERSITY

Appendix A

Reducing the least-squares cubic
equation to commonly known
regression methods under specificassumptions
Let us consider the earlier expression of the LSC slope in Eq. (3.48) (Chapter 3) as written
below,


bVi
ri Vi
Ui
Vi
+
−
Wi2 w(Yαiw(Xi )
i)
 .

b= P
Ui
bVi
Wi2 w(Y+ w(X− brαi Ui i Ui
i)
i)
P

(A.1)

This solution is general, and various other regression methods, such as OLS and ODR, can beobtained only as special cases of this equation. As demonstrated in [13], the sections below
explore through some of the special cases which reduce to some of the commonly known regression methods.

Below is an earlier expression in Eq. (3.35), and here it is noted for eased referencing.
Wi =w(Xi )w(Yi ).
2b w(Yi ) + w(Xi ) − 2bri αi

95

(A.2)

PK SHIGWEDHA

A.1

APPENDIX

Only Yi is subject to errors (i.e., Xi is error-free):In this case, w(Xi ) → ∞, and ri → 0. The expression in Eq. (A.2) then reduces to Wi = w(Yi );and thus Eq. (A.1) solves to,Pw(Yi )Ui Vi
b= P.
w(Yi )Ui2(A.3)

This solution is equivalent to the ussual weighted OLS regression of Y on X. Removing theweight component, Eq. (A.3) then reduces to the expression derived for OLS in Eq. (2.56) as,P
Ui V ib= P 2Ui(A.4)
P(Xi − X̄)(Yi − Ȳ )P=
(Xi − X̄)2

A.2

Only Xi is subject to errors (i.e., Yi is error-free):In this case, w(Yi ) → ∞, and ri → 0. The expression in Eq. (A.2) then reduces to Wi =w(Xi );b2

and thus Eq. (A.1) solves to,Pw(Xi )Vi2.
b= P
w(Xi )Ui Vi(A.5)

This solution is equivalent to the ussual weighted OLS regression of X on Y , and only the X
variable is being weighted. Removing the weight component, it reduces to the OLS equation
below, which is obtained by considering only the measurement errors in the X variable;
P 2
V
b= P iUi V i(A.6)− Ȳ )2P(Yi=P(Xi − X̄)(Yi − Ȳ )A.3

Errors present in both X and Y variables:

If we assume that w(Xi ) = w(Yi ) = c (constant) such that,w(Xi )= 1, and ri = 0, Wi in Eq.
w(Yi )(A.2) becomes,
Wi =cb2 + 1The expression in Eq. (A.1) then reduces to,
2 
P
a

b = P

b2 +1

a
b2 +1

2 

.

(A.7)

Ui
bVi
a + a



Vi

Ui
bVi
a + a



Ui

,

(A.8)


P
Ui Vi + bVi2
=P 2
Ui + bUi Vi
工学院大学

96

KOGAKUIN UNIVERSITYPK SHIGWEDHA

APPENDIX

which gives,b2

X

Ui V i − b

X

Vi2 −

X

 X
Ui2 −
Ui Vi = 0.

(A.9)

Solving for b from Eq. (A.9) using the standard quadratic approach gives,
q P
P 2
P
P
P
Vi2 − Ui2 + 4 ( Ui Vi )2( Vi2 − Ui2 ) +Pb=
2 Ui Vi(A.10)=(Syy − Sxx ) +q

2(Syy − Sxx )2 + 4Sxy

2Sxy

.

Equation (A.10) is identical to Eq. (2.62) for ODR, which is obtained by minimizing the sum of
the squares of the perpendicular distances between the observed data points and the estimatedstraight line. Again it is shown here that this solution is obtained by assuming that the variancesof the X and Y errors are equal, since in cases where the variances are known, the weights (whichwere here are assumed to be equal) are estimated as the inverse of the variances.

工学院大学

97

KOGAKUIN UNIVERSITY

Appendix B

Noise model
The noise model utilised in this study is adopted from previous studies [6, 8, 11, 16]. As per thedescription in [11], given the radioactivity concentration, CT , of the radiotracer in the ROI, if
c is the count measured at mid-time t of the frame with ”frame length” ∆t, then we have that,c=

CT · d · ∆t,awhere a is a scaling factor, and d is a decay factor given by, d =

(B.1)

1 t/T
, where T is the half-life2of a radionuclide. The error in CT can be expressed as,∆CT =where, ∆c =√a · ∆c,d · ∆t(B.2)

c, is the error in c. So, the simulated noisy tTACs are then given by
C¯T = CT (1 + cv · G(0, 1)),

(B.3)

r
∆CT
a
where, cv ==, is the variation coefficient, and G(0, 1) is a random numberCT
CT · d · ∆tof a Gaussian distribution with mean 0, and standard deviation 1.

98

Appendix C

Contrast evaluation of simulationdata
For contrast comparison between LSC–PCA and OLS–PCA in simulation data for 11 C-PiB, a
set of hypothetical white matter regions were simulated —in addition to the gray matter regions
simulated in Section 4.1.1. For 11 C-PiB Aβ-negative images, white matter regions can have highretention of the radiotracer due to nonspecific retention [25, 56], and slow accumulation andslower clearance [25]. Therefore, the white matter regions were simulated with the values of K1(0.190 mL cm−1 ) and k2 (0.100 min−1 ) lower than those of the gray matters, and DVR values
high than those of the gray matters. Thus, to simulate white matter with high DVR values
than the garay matters, the 11 unit vector, [0.4 : 0.055 : 0.95], was added to the 11 unit vectorof the gray matter DVRs. This gave the lowest white matter DVR a value of 1.6727 (the lowest
for gray matter is 1.2727), and the highest of 4.95 (the highest for gray matter is 4). The ratioof these gray and white matter DVR values are within the ranges reported in a previous study[57]. The scan period and number of frames for 11 C-PiB white matter simulations were set tobe similar to those of the gray matter, and thus similar to those of clinical 11 C-PiB PET images
used in this study.

Similar with gray matter simulations, a noise-free tTAC was formed for each white matter DVR
value. Statistical noise was added to the noise-free tTACs to form noisy tTACs. For eachnoise-free tTAC, 1024 noisy tTACs were simulated as a slice of 32 × 32 voxels. Four methods,
LSC–PCA, OLS–PCA, LSC, and conventional OLS-based LGA, were then used to estimate the
DVR values from the noisy tTACs, and the results were compared.

99

PK SHIGWEDHA

APPENDIX

The simulated gray and white matter data for 11 C-PiB were denoised together by PCA. For
each data set (1024 simulated tTACS for each of the 11 true DVRs), DVR estimation from thenoisy tTACs were performed thrice at 9 percentage variance level of PCA, 91%, 92%, 93%, 94%,
95%, 96%, 97%, 98%, 99% (i.e., 27 times) for for LSC–PCA and OLS–PCA methods. This wasto assess the performance of the two methods over a wide range of percentage variance, sincein a real clinical practice the actual percentage variance to be retained in PCA denoising is notknown. The three estimates at each variance level were differed by a slight change of < 5%in the magnitude of noise level. The magnitude of the noise level refers to the percentage ratio of the means of “standard deviation in the noisy tTACs to the mean of the noise-free tTAC”.

To have equal set of data for all methods, the LSC and OLS estimates were performed 27 times
at noise levels corresponding to those of the PCA-based methods. The mean and standard
deviations were then calculated for 27 estimates for each 1024 simulated tTACs for each true
DVR value.

Figure C.1 compares the results of contrast between the simulated gray and white matter regions. The comparisons include all calculated DVR values (three times at 9 percentage variancelevels for LSC–PCA and OLS–PCA, and at 27 different noise levels for LSC and OLS). Statistical comparisons of these contrast are shown in table C.1, using the Friedman test [58] andPost-Hoc multiple comparison. A statistical significant difference is observed between LSC–
PCA and all other methods. On the other hand, no statistical significant difference is observed
between OLS–PCA and LSC.

The significant p value observed between LSC–PCA and OLS–PCA in simulation results could
be an indication that the same can be obtained for clinical data. This can be investigated witha larger cohort of clinical 11 C-PiB brain data.

工学院大学

100

KOGAKUIN UNIVERSITYPK SHIGWEDHA

References'$

LSC-PCAOLS-PCALSC
OLSContrast

1.5

1

0.5

0
1.25 1.525 1.8 2.0752.35 2.625 2.9 3.1753.45 3.725 4&%Figure C.1: Box plots of the DVR contrast between the hypothetical simulated gray and whitematter regions, calculated for all DVR values estimated at different percentage variance levels
for LSC–PCA and OLS–PCA, and at 27 different noise levels for LSC and OLS. The circles inthe box plots denote the means, whereas the horizontal lines denote the medians.

Table C.1: Friedman test and Post-Hoc multiple comparisons of the p-values of the contrast
values of simulated 11 C-PiB data plotted in figure C.1. Tabulated numerical values correspondto the p-values between the two methods in the corresponding rows and columns.
Friedman Test (p value < 0.001)
Post-Hoc Multiple Comparisons p-values

LSC–PCAOLS-PCA

LSC

OLS

0.041

0.022

< 0.001

0.804< 0.009

OLS–PCA
LSC

工学院大学

0.018

101

KOGAKUIN UNIVERSITY

Bibliography[1] D. L. Bailey, M. N. Maisey, D. W. Townsend, and P. E. Valk, Positron Emission Tomography: Basic Sciences. Springer, 2005.
[2] J. Logan, J. S. Fowler, N. D. Volkow, A. P. Wolf, S. L. Dewey, D. J. Schlyer, R. R. MacGregor, R. Hitzemann, B. Bendriem, S. J. Gatley, et al., “Graphical analysis of reversible radioligand binding from time—activity measurements applied to [N-11 C-methyl]-(–)-Cocaine
PET studies in human subjects,” Journal of Cerebral Blood Flow & Metabolism, vol. 10,
no. 5, pp. 740–747, 1990.
[3] J. Logan, J. S. Fowler, N. D. Volkow, G.-J. Wang, Y.-S. Ding, and D. L. Alexoff, “Distribution volume ratios without blood sampling from graphical analysis of pet data,” Journalof Cerebral Blood Flow & Metabolism, vol. 16, no. 5, pp. 834–840, 1996.
[4] M. Slifstein and M. Laruelle, “Effects of statistical noise on graphic analysis of pet neuroreceptor studies,” Journal of Nuclear Medicine, vol. 41, no. 12, pp. 2083–2088, 2000.
[5] J. Logan, “A review of graphical methods for tracer studies and strategies to reduce bias,”
Nuclear Medicine and Biology, vol. 30, no. 8, pp. 833–844, 2003.
[6] Y. Kimura, M. Naganawa, M. Shidahara, Y. Ikoma, and H. Watabe, “PET kinetic analysis—pitfalls and a solution for the logan plot,” Annals of Nuclear Medicine, vol. 21, no. 1,
pp. 1–8, 2007.
[7] A. Abi-Dargham, D. Martinez, O. Mawlawi, N. Simpson, D.-R. Hwang, M. Slifstein, S. Anjilvel, J. Pidcock, N.-N. Guo, I. Lombardo, et al., “Measurement of striatal and extrastriatal
dopamine d1 receptor binding potential with [11 C] NNC 112 in humans: validation andreproducibility,” Journal of Cerebral Blood Flow & Metabolism, vol. 20, no. 2, pp. 225–243,
2000.

102

PK SHIGWEDHA

References[8] J. Logan, J. S. Fowler, N. D. Volkow, Y. S. Ding, G.-J. Wang, and D. L. Alexoff, “A strategy
for removing the bias in the graphical analysis method,” Journal of Cerebral Blood Flow
& Metabolism, vol. 21, no. 3, pp. 307–320, 2001.
[9] M. Ichise, H. Toyama, R. B. Innis, and R. E. Carson, “Strategies to improve neuroreceptor
parameter estimation by linear regression analysis,” Journal of Cerebral Blood Flow &Metabolism, vol. 22, no. 10, pp. 1271–1281, 2002.
[10] M. Ichise, J.-S. Liow, J.-Q. Lu, A. Takano, K. Model, H. Toyama, T. Suhara, K. Suzuki,
R. B. Innis, and R. E. Carson, “Linearized reference tissue parametric imaging methods:
application to [11 C] DASB positron emission tomography studies of the serotonin transporter in human brain,” Journal of Cerebral Blood Flow & Metabolism, vol. 23, no. 9,
pp. 1096–1112, 2003.
[11] J. Varga and Z. Szabo, “Modified regression model for the logan plot,” Journal of Cerebral
Blood Flow & Metabolism, vol. 22, no. 2, pp. 240–244, 2002.
[12] D. York, N. M. Evensen, M. L. Martı́nez, and J. De Basabe Delgado, “Unified equations
for the slope, intercept, and standard errors of the best straight line,” American Journalof Physics, vol. 72, no. 3, pp. 367–375, 2004.
[13] D. York, “Least-squares fitting of a straight line,” Canadian Journal of Physics, vol. 44,
no. 5, pp. 1079–1086, 1966.
[14] D. York, “Least squares fitting of a straight line with correlated errors,” Earth and Planetary Science Letters, vol. 5, pp. 320–324, 1968.
[15] I. T. Jolliffe and J. Cadima, “Principal component analysis: a review and recent developments,” Philosophical Transactions of the Royal Society A: Mathematical, Physical andEngineering Sciences, vol. 374, no. 2065, p. 20150202, 2016.
[16] A. Joshi, J. A. Fessler, and R. A. Koeppe, “Improving pet receptor binding estimatesfrom logan plots using principal component analysis,” Journal of Cerebral Blood Flow &Metabolism, vol. 28, no. 4, pp. 852–865, 2008.
[17] Y. Kimura, M. Senda, and N. M. Alpert, “Fast formation of statistically reliable FDG
parametric images based on clustering and principal components,” Physics in Medicine &Biology, vol. 47, no. 3, p. 455, 2002.
工学院大学

103

KOGAKUIN UNIVERSITYPK SHIGWEDHA

References[18] J. Axelsson and J. Sörensen, “The 2d hotelling filter-a quantitative noise-reducing principalcomponent filter for dynamic PET data, with applications in patient dose reduction,” BMC
Medical Physics, vol. 13, no. 1, pp. 1–17, 2013.
[19] L. C. Parra, S. Haufe, and J. P. Dmochowski, “Correlated components analysis-extracting
reliable dimensions in multivariate data,” arXiv preprint arXiv:1801.08881, 2018.
[20] P. K. Shigwedha, T. Yamada, K. Hanaoka, K. Ishii, Y. Kimura, and Y. Fukuoka, “Astrategy to account for noise in the x-variable to reduce underestimation in logan graphicalanalysis for quantifying receptor density in positron emission tomography,” BMC Medical
Imaging, vol. 20, no. 1, pp. 1–8, 2020.
[21] H. Watabe, Y. Ikoma, Y. Kimura, M. Naganawa, and M. Shidahara, “PET kinetic analysis—compartmental model,” Annals of Nuclear Medicine, vol. 20, no. 9, pp. 583–588,
2006.
[22] A. Joshi, Improved quantitative methods for multiple neuropharmacological non-invasive
brain PET studies. PhD thesis, University of Michigan, 2008.
[23] E. Miele, G. P. Spinelli, F. Tomao, A. Zullo, F. De Marinis, G. Pasciuti, L. Rossi, F. Zoratto,
and S. Tomao, “Positron emission tomography (PET) radiotracers in oncology–utility of
18 F-fluoro-deoxy-glucose (FDG)-PET in the management of patients with non-small-celllung cancer (NSCLC),” Journal of Experimental & Clinical Cancer Research, vol. 27, no. 1,
pp. 1–10, 2008.
[24] S. Yu, “Review of 18 F FDG synthesis and quality control,” Biomedical Imaging and Intervention Journal, vol. 2, no. 4, 2006.
[25] W. E. Klunk, H. Engler, A. Nordberg, Y. Wang, G. Blomqvist, D. P. Holt, M. Bergström,
I. Savitcheva, G.-F. Huang, S. Estrada, et al., “Imaging brain amyloid in alzheimer’s disease
with pittsburgh compound-b,” Annals of Neurology: Official Journal of the American
Neurological Association and the Child Neurology Society, vol. 55, no. 3, pp. 306–319,
2004.
[26] J. C. Price, W. E. Klunk, B. J. Lopresti, X. Lu, J. A. Hoge, S. K. Ziolko, D. P. Holt,
C. C. Meltzer, S. T. DeKosky, and C. A. Mathis, “Kinetic modeling of amyloid binding inhumans using PET imaging and pittsburgh compound-b,” Journal of Cerebral Blood Flow
& Metabolism, vol. 25, no. 11, pp. 1528–1547, 2005.
工学院大学

104

KOGAKUIN UNIVERSITYPK SHIGWEDHA

References[27] J. J. Vaquero and P. Kinahan, “Positron emission tomography: current challenges and opportunities for technological advances in clinical and preclinical imaging systems,” Annual
Review of Biomedical Engineering, vol. 17, pp. 385–414, 2015.
[28] A. K. Shukla and U. Kumar, “Positron emission tomography: An overview,” Journal of
Medical Physics/Association of Medical Physicists of India, vol. 31, no. 1, p. 13, 2006.
[29] E. D. Morris, C. J. Endres, K. C. Schmidt, B. T. Christian, R. F. Muzic, and R. E.
Fisher, “Kinetic modeling in positron emission tomography,” Emission Tomography: The
Fundamentals of PET and SPECT. Academic, San Diego, 2004.
[30] R. Koeppe, V. Holthoff, K. Frey, M. Kilbourn, and D. Kuhl, “Compartmental analysis of[11 C] flumazenil kinetics for the estimation of ligand transport rate and receptor distribution using positron emission tomography,” Journal of Cerebral Blood Flow & Metabolism,
vol. 11, no. 5, pp. 735–744, 1991.
[31] Y. Ikoma, H. Watabe, M. Shidahara, M. Naganawa, and Y. Kimura, “PET kinetic analysis: error consideration of quantitative analysis in dynamic studies,” Annals of Nuclear
Medicine, vol. 22, no. 1, pp. 1–11, 2008.
[32] M. Ichise, J. R. Ballinger, H. Golan, D. Vines, A. Luong, S. Tsai, and H. F. Kung, “Noninvasive quantification of dopamine d2 receptors with iodine-123IBF [SPECT,” Journal of
Nuclear Medicine, vol. 37, no. 3, pp. 513–520, 1996.
[33] M. Yaqub, N. Tolboom, R. Boellaard, B. N. van Berckel, E. W. van Tilburg, G. Luurtsema, P. Scheltens, and A. A. Lammertsma, “Simplified parametric methods for [11 C] pib
studies,” Neuroimage, vol. 42, no. 1, pp. 76–86, 2008.
[34] R. B. Innis, V. J. Cunningham, J. Delforge, M. Fujita, A. Gjedde, R. N. Gunn, J. Holden,
S. Houle, S.-C. Huang, M. Ichise, et al., “Consensus nomenclature for in vivo imaging ofreversibly binding radioligands,” Journal of Cerebral Blood Flow & Metabolism, vol. 27,
no. 9, pp. 1533–1539, 2007.
[35] M. A. Mintun, M. E. Raichle, M. R. Kilbourn, G. F. Wooten, and M. J. Welch, “A quantitative model for the in vivo assessment of drug binding sites with positron emission tomography,” Annals of Neurology: Official Journal of the American Neurological Associationand the Child Neurology Society, vol. 15, no. 3, pp. 217–227, 1984.
工学院大学

105

KOGAKUIN UNIVERSITYPK SHIGWEDHA

References[36] L. Farde, L. Eriksson, G. Blomquist, and C. Halldin, “Kinetic analysis of central [11 C]
raclopride binding to d2-dopamine receptors studied by PET—a comparison to the equilibrium analysis,” Journal of Cerebral Blood Flow & Metabolism, vol. 9, no. 5, pp. 696–708,
1989.
[37] Y. Zhou, W. Ye, J. R. Brašić, A. H. Crabb, J. Hilton, and D. F. Wong, “A consistent andefficient graphical analysis method to improve the quantification of reversible tracer bindingin radioligand receptor dynamic pet studies,” Neuroimage, vol. 44, no. 3, pp. 661–670, 2009.
[38] D. Feng, S.-C. Huang, Z. Wang, and D. Ho, “An unbiased parametric imaging algorithm
for nonuniformly sampled biomedical system parameter estimation,” IEEE Transactions
on Medical Imaging, vol. 15, no. 4, pp. 512–518, 1996.
[39] X. Dai, Z. Chen, and J. Tian, “Performance evaluation of kinetic parameter estimation
methods in dynamic FDG-PET studies,” Nuclear Medicine Communications, vol. 32, no. 1,
pp. 4–16, 2011.
[40] R. Wehr and S. R. Saleska, “The long-solved problem of the best-fit straight line: Application to isotopic mixing lines,” Biogeosciences, vol. 14, no. 1, pp. 17–29, 2017.
[41] J. Fox, Applied regression analysis and generalized linear models. Sage Publications, 3 ed.,2015.
[42] J. M. Wooldridge, Introductory Econometrics: A Modern Approach. Cengage Learning,
6 ed., 2015.
[43] W. M. Stock J, Introduction to Econometrics. Pearson Addison Wesley, 2008.
[44] P. T. Boggs and J. E. Rogers, “Orthogonal distance regression,” Contemporary Mathematics, vol. 112, pp. 183–194, 1990.
[45] R. Carroll and D. Ruppert, “The use and misuse of orthogonal regression in linear errorsin-variables models,” The American Statistician, vol. 50, no. 1, pp. 1–6, 1996.
[46] R. Deakin and M. Hunter, “Fitting a line of best fit to correlated data of varying precision.” http://www.mygeodesy.id.au/documents/Fitting\%20a\%20straight\%20line.pdf, 2017. Accessed: 2019-05-27.
[47] W. E. Deming, Statistical Adjustment of Data. Wiley, 1943.
工学院大学

106

KOGAKUIN UNIVERSITYPK SHIGWEDHA

References[48] N. Shammas, “Hp 39gii regression: Part iv least-squares relative error regression.” http://www.namirshammas.com/HP39gII/HP39gII\_part4.pdf, 2013. Accessed: 2018-12-26.
[49] C. Tofallis, “Least squares percentage regression,” Journal of Modern Applied Statistical
Methods, vol. 7, no. 2, pp. 526–534, 2008.
[50] J. O. Rawlings, S. G. Pantula, and D. A. Dickey, Applied regression analysis: a researchtool. New York: Springer, 2 ed., 2001.
[51] A. Cottrell, “Regression analysis: basic concepts.” http://users.wfu.edu/cottrell/
ecn215/regress.pdf. Accessed: 2018-12-26.
[52] C. J. Endres, B. Bencherif, J. Hilton, I. Madar, and J. J. Frost, “Quantification of brainµ-opioid receptors with [11 C] carfentanil: reference-tissue methods,” Nuclear Medicine andBiology, vol. 30, no. 2, pp. 177–186, 2003.
[53] I. Jolliffe, “Principal component analysis, 2nd edn. series: Springer series in statistics, vol.
XXIX, 487,” Springer, NY, p. 28, 2002.
[54] H. Abdi and L. J. Williams, “Principal component analysis,” Wiley Interdisciplinary Reviews: Computational Statistics, vol. 2, no. 4, pp. 433–459, 2010.
[55] T. Yamada, Y. Kimura, K. Fujii, S. Watanabe, T. Nagaoka, M. Nemoto, K. Hanaoka,
H. Kaida, C. Hosokawa, and K. Ishii, “Noise reduction algorithm for amyloid imaging
preserving image resolution–quantitative evaluation using clinical images,” Human Amyloid
Imaging Conference, 2018.
[56] R. A. Sperling, P. S. LaViolette, K. O’Keefe, J. O’Brien, D. M. Rentz, M. Pihlajamaki,
G. Marshall, B. T. Hyman, D. J. Selkoe, T. Hedden, et al., “Amyloid deposition is associated with impaired default network function in older persons without dementia,” Neuron,vol. 63, no. 2, pp. 178–188, 2009.
[57] C. Hosokawa, K. Ishii, Y. Kimura, T. Hyodo, M. Hosono, K. Sakaguchi, K. Usami, K. Shimamoto, Y. Yamazoe, and T. Murakami, “Performance of 11c-pittsburgh compound b petbinding potential images in the detection of amyloid deposits on equivocal static images,”
Journal of Nuclear Medicine, vol. 56, no. 12, pp. 1910–1915, 2015.

工学院大学

107

KOGAKUIN UNIVERSITYPK SHIGWEDHA

References[58] M. Friedman, “The use of ranks to avoid the assumption of normality implicit in the analysis
of variance,” Journal of the american statistical association, vol. 32, no. 200, pp. 675–701,
1937.

工学院大学

108

KOGAKUIN UNIVERSITY