Jukuri, open repository of the Natural Resources Institute Finland (Luke) 
   
 
   
All material supplied via Jukuri is protected by copyright and other intellectual property rights. Duplication 
or sale, in electronic or print form, of any part of the repository collections is prohibited. Making electronic 
or print copies of the material is permitted only for your own personal use or for educational purposes.  For 
other purposes, this article may be used in accordance with the publisher’s terms. There may be 
differences between this version and the publisher’s version. You are advised to cite the publisher’s 
version. 
 
This is an electronic reprint of the original article.  
This reprint may differ from the original in pagination and typographic detail. 
 
Author(s): Mikko Kuronen, Mari Myllymäki, Adam Loavenbruck and Aila Särkkä 
Title: Point process models for sweat gland activation observed with noise 
Year: 2021 
Version: Published version 
Copyright:   The Author(s) 2021 
Rights: CC BY 4.0 
Rights url: http://creativecommons.org/licenses/by/4.0/ 
 
Please cite the original version: 
Kuronen, M, Myllymäki, M, Loavenbruck, A, Särkkä, A. Point process models for sweat gland 
activation observed with noise. Statistics in Medicine. 2021; 1– 18. 
https://doi.org/10.1002/sim.8891 
Received: 27 April 2020 Revised: 7 November 2020 Accepted: 6 January 2021
DOI: 10.1002/sim.8891
RE S EARCH ART I C L E
Point process models for sweat gland activation observed
with noise
Mikko Kuronen1 Mari Myllymäki1 Adam Loavenbruck2 Aila Särkkä3
1Natural Resources Institute Finland
(LUKE), Helsinki, Finland
2Department of Neurology, Kennedy
Laboratory, University of Minnesota,
Minneapolis, Minnesota
3Department of Mathematical Sciences,
Chalmers University of Technology and
the University of Gothenburg,
Gothenburg, Sweden
Correspondence
Mikko Kuronen, Natural Resources
Institute Finland (LUKE), Helsinki,
Finland.
Email: mikko.p.kuronen@gmail.com
Funding information
Academy of Finland, Grant/Award
Numbers: 295100, 306875, 32711; Swedish
Foundation for Strategic Research,
Grant/Award Number: SSF AM13-0066;
Swedish Research Council, Grant/Award
Number: VR 2013-5212
The aim of this article is to construct spatial models for the activation of sweat
glands for healthy subjects and subjects suffering from peripheral neuropathy
by using videos of sweating recorded from the subjects. The sweat patterns are
regarded as realizations of spatial point processes and two point process mod-
els for the sweat gland activation and two methods for inference are proposed.
Several image analysis steps are needed to extract the point patterns from the
videos and some incorrectly identified sweat gland locations may be present in
the data. To take into account the errors, we either include an error term in the
point process model or use an estimation procedure that is robust with respect
to the errors.
KEYWORD S
Bayesian inference, independent thinning, peripheral neuropathy, point pattern, sequential point
process, soft-core inhibition
1 INTRODUCTION
Sweating is critical in human evolution in maintaining ability to thermoregulate in a wide range of climates and activity
levels. Neurologic control, headquartered in the hypothalamis, is therefore tightly regulated and concordantly disruptive
in pathologic states such as peripheral neuropathy. Assessment of sudomotor (sweat) function has long been used in
clinical and research settings for detection of neurologic dysfunction.1,2 Minor’s starch iodine test was originally described
in 1928,3 where tincture of iodine was applied to the skin, air dried, and then powdered with corn starch. Sweating is
stimulated with increasing room temperature or use of pilocarpine. As sweat flows from pores, iodine is diluted and the
solution absorbed by the starch powder, turning dark black from yellow. Normally the entire skin surface should be able
to sweat in response to sufficient stimuli, and absence of sweating in an area of the body is indicative of loss of neurologic
function in that area.
Peripheral neuropathy is a disease state of peripheral nerves, the segment of the nervous system which extends from
the brain and spinal cord to various targets in the body, such as muscles, sensory receptors, and autonomically controlled
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the
original work is properly cited.
© 2021 The Authors. Statistics in Medicine published by John Wiley & Sons Ltd.
Statistics in Medicine. 2021;1–18. wileyonlinelibrary.com/journal/sim 1
2 KURONEN et al.
organs. Peripheral neuropathy occurs in etiologically diverse conditions which cause damage or dysgenesis of peripheral
nerves. The most common causes include diabetes, toxicity such as in chemotherapy and excessive alcohol consumption,
and vitamin deficiencies.4 The resulting nerve damage causes various combinations ofmuscleweakness, pain, numbness,
and autonomic dysfunction.
Autonomic nerves are often the earliest to be affected in peripheral neuropathy,5-7 and penetrate all parts of the body,
including digestive tract, liver, kidneys, bladder, genitals, lungs, pupils, heart, and skin. Skin is the largest organ in the
body, and contains a vast network of the distal ends of sensory and autonomic nerves over the entire body surface.
These distal ends of nerves are especially susceptible to systemic disease. Because sweating is neutrally controlled and
modulated, and can be measured at the skin surface, it can be used to reflect alterations in the underlying nerves.
Currently, there are several tests used in clinical practice to evaluate sudomotor function.8-10 Thermoregulatory sweat
testing11 measures percentage of body surface area sweating elicited by a heated, humidified sauna. Sweat imprints and
silastic molds12,13 measure the density and distribution of activated sweat glands in an area of skin. Quantitative sudomo-
tor axon reflex testing (QSART) is likely the most widely used autonomic test of sweating,14,15 comparing against robust
normative databases the total volume of sweat produced by 1 cm2 areas of skin at standardized sites.
The sensitive sweat test (SST) enhanced Minor’s starch iodine test with closeup time lapse imaging, and software
analysis.4,16-18 The critical feature of the test is a rigid, transparent video screen which limits sweat to an essentially
two-dimensional space. As sweat exits ducts, it dilutes the iodine painted on the skin onto starch coated plastic film. The
imaged result is a field of sharply demarcated, dark sweat spots on a white background, expanding centripetally around
the opening of each duct (Figure 1). The area of each spot is therefore a measurement of the volume of sweat produced
by each gland. Sub-nanoliter volumes of liquid were measured by pipette and shown to create reproducible sweat spot
areas. Similarly, tracking the increase in sweat spots’ areas between timelapse frames measures the rate of sweating from
each duct at the nanoliter level. Of note, blackened areas of film do not return to white during the test—sweat spots can
only enlarge. The videos therefore provide several measurable physiologic data points, including coordinates and rela-
tive locations of all sweat spots, the second by second volumes of sweat (nanoliters) and flow rate of each sweat gland
(nL/minute), total number of activated sweat glands, density of activated sweat glands (glands/cm2), total sweat volume
(nL), and total sweat rate (nL/minute).
Using the dynamic sweat test, a significant reduction of sweating was observed in diabetic subjects in the distal leg
but not in forearm.16 The study included measurements taken from the forearm of 14 diabetic subjects and 14 age- and
sex-matched healthy controls and from distal leg of seven diabetes patients and seven controls. In a larger study,4 178
healthy controls and 20 neuropathy subjects were tested, most of them at the hand, thigh, calf, and foot, some only at calf
and foot, and it was concluded that neuropathy subjects had lower sweat rates per sweat gland, lower total amount of
sweat, and lower sweat gland density. It was also observed visually that the sweat patterns of the diabetic subjects were less
regular than the healthy patterns.16 To quantify this visual observation, the spatial sweat patterns that the videos provide
should be investigated more carefully. Spatial analysis can provide more detailed information on the sweat patterns and
tools for revealing additional abnormalities that may appear when the sweat rate and the total amount of sweat are still
within normal ranges. However, up to now, no spatial analysis of the sweat patterns to quantify this observation has been
performed.
In this article, the locations of sweat spots or glands are regarded as realizations of spatial point processes. Our main
emphasis is to develop suitable methodology for analyzing the spatial structure of the sweat gland patterns and activation
extracted from video sequences. As the data are non-standard in the point process literature, some special treatment is
needed.
To extract the coordinates of the individual sweat spots, that is, the points of the point patterns, from the videos (see
Figure 1), several image analysis steps are needed. As a non-standard step, we apply an algorithmbased on the detection of
a change point to the sequence of gray scale values in each pixel. This pixel-by-pixel approach suits to the video sequences,
where the sweat does not dry once it has appeared, much better than going through the videos frame by frame, because
the sweat gland locations are best detected at times where sweat first appears. However, even though we perform careful
analysis of the videos, there are some spots that are incorrectly recorded as two spots due to, for example, wrinkles in the
skin. It is not straightforward to remove these errors automatically and doing it manually can be very time consuming.
Therefore, they need to be taken into account in the analysis of the point patterns.
Some studies of point patterns observed with errors or noise can be found in the literature. For example, in the area of
minefield detection, one first detects a minefield and then classifies each observed point in the minefield either as mine
or as noise. The observation window is typically divided into two parts, the minefield as a region with a higher intensity
containing both mines and noise and a low intensity area containing only noise.19 The points can then be classified in a
KURONEN et al. 3
F I GURE 1 A sequence of snapshots (1 second (top left), 15 seconds (top right), 30 seconds (bottom left), and 60 seconds (bottom
right)) of one control subject with extracted gland locations (+)
Bayesian set-up where a posterior probability for each point being amine is derived.20,21 In a similar Bayesian framework,
classification of points of a superposition of a Strauss process and Poisson noise has been considered.22,23 A likelihood of
an imperfect observation given the true point process, where the imperfect observation can be due to random thinning,
displacement, censoring of the displaced points, or superposition of additional points24 and Bayesian analysis for similar
data25 can also be found in the literature. Furthermore, a Bayesian framework for estimating the intensity of a noisy point
process, where the noise is either due to points within the sampling window but regarded as being outside and/or points
outside the window which were incorrectly regarded as points inside the window, is available26 as well as descriptive
statistics for noisy spatial point processes, where the noise is perturbation of points.27
Here, we suggest two different ways to model the activation of sweat glands and to take noise into account in the
analysis, by either including an error term in the model or using an estimation procedure that is robust with respect to
the errors. We pay special attention to incorrectly recorded close pairs of points since they can cause problems for the
analysis of regular point patterns such as our data.
We propose two models for the sweat gland patterns that are different in their nature. In the first model, the activa-
tion of individual sweat glands is modeled by a sequential point process, where sweat spots appear conditionally on the
pattern so far. The other model is more physiologically motivated, a generative point process model, where the activated
sweat gland pattern is modeled as a thinning of the underlying true (unobserved) sweat gland pattern which is modeled
first.While the likelihood of the sequential model is tractable, it has been considered computationally costly to evaluate.28
Here, we propose an efficient way to perform traditional likelihood-based inference for a certain type of sequential mod-
els, which makes also likelihood-based Bayesian inference feasible. The likelihood of the generative model is not easily
tractable and, therefore, we employ approximate Bayesian computation (ABC) to estimate the model parameters. When
some noise points are present, the sequential model is replaced by a mixture model having the sequential point process
and an error point process as its components. In the generative model, the summary statistics in the ABC approach are
chosen such that they are robust with respect to the errors.
4 KURONEN et al.
The rest of this article is organized as follows. We first describe the extraction of the points from the videos and the
preliminary analysis of the data in Section 2. Then, we present the sequential and generative models together with a
description of the inferencemethods in Sections 3 and 4, respectively. Further, themethodology is illustrated by analyzing
a set of video sequences taken from the right foot of 15 subjects, either controls or subjects with suspected or diagnosed
neuropathy. The models are fitted separately to each subject. Section 5 is left for further discussion. We carry out all
computations in Julia29 while we use R30 mainly for plotting.
2 DATA AND PRELIMINARY DATA ANALYSIS
2.1 Description of data
The data have been collected by Dr. Kennedy’s group at the University of Minnesota by using the dynamic sweat test they
have presented.16 First, sweating is stimulated by placing a patch with pilocarpine gels on the test site, foot or calf. Then,
the test site is dried and painted with iodine solution. Finally, a camera is placed on the skin and a video is recorded at the
rate of 1 frame/sec for 60 seconds. The size of the frame was 2592× 1944 pixels corresponding to 17.5× 13 mm2. Videos
were taken from the feet and/or calves from 121 healthy controls without known neuropathy or known risk factors for
neuropathy, as well as 72 subjects who had reported having symptoms of neuropathy, 20 of whom had well characterized
neuropathy (diagnosis based on neurologic examination and nerve conduction studies). Therefore, the subjects were
divided into three groups: controls, subjects with suspected neuropathy (MNA), and subjects with diagnosed neuropathy
(MNA Diagnosed).
In this study, we have access to five videos from the right foot from each of the three groups. Based on earlier studies,
it was clear that the number of activated sweat glands is an important predictor for the condition, controls having higher
density than subjects in the neuropathy groups. The five videos were selected based on the point density of the pattern
so that different groups have overlapping densities. A sequence of snapshots (1 second, 15 seconds, 30 seconds, and 60
seconds) of one control subject is shown in Figure 1. Here, we study the patterns of activated sweat glands at the end
of videos as realizations of spatial point processes. The complete video is needed to extract the gland locations, because
these can be obtained most precisely at their first occurrence (see Section 2.2).
2.2 Video processing with change point detection
Extracting the locations and sizes of the sweat spots required several video processing steps: transforming the video
into sweat/not sweat binary video, splitting the sweat part of the video into the sweat produced by individual sweat
glands, and finally extracting the point pattern of gland locations. As is commonly done in data pre-processing, some
trial and error was necessary before a satisfactory result was obtained. In the following, we describe the final choices in
more detail.
The first step consisted of a background correction, finding change points, and finding and applying a threshold to the
change points. The backgroundwas first estimated by kernel smoothing using aGaussian kernelwith 𝜎 = 100 pixels to the
first frame of the video. The smoothing bandwidth of 100 pixels, which is a bitmore than theminimum interpoint distance
of sweat spots, was large enough to remove the sweat spots, but small enough to keep the background fluctuations. Since
the first frame had only small amount of sweat, the resulting imagemostlymimicked the lighting conditions. For example,
the corners of the frame were darker than the middle. Next, the grayscale values gt of each pixel at frame t were divided
by the estimated lightning intensity l of the pixel and the time series of these scaled grayscale values were considered to
find the pixels that belong to the wet area. More precisely, a time series was constructed for each pixel as follows: Let
x−2 = x−1 = x0 = 1 and xt = gt/l for t= 1, … ,T. The change point of the time series x−2, x−1, x0, x1, … , xT was defined as
the integer value t≥ 1 that minimizes f (t) = s2−2,t + s2t+1,T , where
s2j,k =
1
k − j + 1
k∑
i=j
x2i −
(
1
k − j + 1
k∑
i=j
xi
)2
is the sample variance of xj, … , xk.
KURONEN et al. 5
F I GURE 2 Time series for four
pixels with estimated jump locations
(frames) marked by dashed lines [Color
figure can be viewed at
wileyonlinelibrary.com]
0.6
0.7
0.8
0.9
1.0
0 20 4 60 0
t
x t
1 2 3 4
Testing for a change point is a well known problem in statistics.31-34 Here, however, the problem is not of purely
statistical nature. There appeared to be some jumps possibly due to the lightning conditions that were not real jumps but
clear enough to be detected by a statistical change point test. In Figure 2, the series number 3 (blue curve) is an example
of such a change point. We used a similar principle as in the statistical change point tests to locate the most likely change
point. Instead of formally testing for the change point, we only accepted those change points where the change was large
enough.
The time series and estimated change points for four pixels are shown in Figure 2. Since each pixel, even the ones that
do not belong to the wet area, obtained a change point, thresholding on the difference of sample means was used to filter
out small changes. A per video threshold was found by trial and error evaluating the point patterns that resulted from
the whole video processing visually. By looking at the video it was quite easy to see the emerging sweat glands and how
well the detected points matched them. In Figure 2, the two largest jumps, 1 and 2, passed the threshold. The third jump,
although it clearly is a jump, did not pass the threshold. The resulting binary video frames were post processed with a
morphological closing to fill in some small gaps.
The sweat area in the first framewas segmented into the sweat produced by individual glands. Startingwith the second
frame, each new sweat pixel was assigned to the closest spot in the previous frame. The distance was measured as the
shortest path through the new sweat area. Several filtering steps were applied in various stages of the process to account
for pixels that belonged to spots that were too small to be sweat.
Finally, we extracted a point pattern with coordinates for each gland. To obtain an ordered point pattern we used the
frame of the first appearance, and for those spots that arrived at the same frame, we used spot size as a surrogate for the
time, where larger ones were assumed to have appeared earlier. An example of extracted point patterns in the video can
be seen in Figure 1. Figure 3 shows the final point patterns of all subjects.
2.3 Spatial summary functions
To analyze the spatial structure of the activated sweat gland patterns, we used two different commonly used spatial
summary functions, the pair correlation function g and the empty space function F. If the underlying point process is
stationary and isotropic, these summary functions are functions of distance r only.
The pair correlation function g describes the (pairwise) interaction between the points.35 It gives a scaled measure
that describes how likely two points are to occur at distance r from each other. For a completely spatially random point
pattern, g(r)= 1 for all r. Values g(r)< 1 indicate inhibition and values g(r)> 1 attraction at distance r. Thus, the pair
correlation function can recognize clustering and regularity at different spatial scales, and is especially useful in describ-
ing strongly inhibitive point patterns. To estimate the pair correlation function, we used a traditional kernel smoothing
method with Epanechnikov kernel and the recommended bandwidth 0.15∕
√
?̂?, where ?̂? is the intensity estimated from
the point pattern.36 Formally,
ĝ(r) = 1
2𝜋r?̂?2
n∑
i=1
n∑
j=1;j≠i
w(xi, xj)k(r − ||xi − xj||),
6 KURONEN et al.
MNA Diagnosed, 23 MNA Diagnosed, 36 MNA Diagnosed, 42 MNA Diagnosed, 50 MNA Diagnosed, 73
MNA, 10 MNA, 20 MNA, 40 MNA, 61 MNA, 71
Controls, 96 Controls, 97 Controls, 149 Controls, 203 Controls, 205
0 1000 2000 0 1000 2000 0 1000 2000 0 1000 2000 0 1000 2000
0
500
1000
1500
2000
0
500
1000
1500
2000
0
500
1000
1500
2000
x [px]
y
[p
x]
F IGURE 3 Point patterns extracted from the videos for each subject of the three groups (group and subject number shown on label)
where n is the observed number of points in the observation window, k is the chosen kernel function, ||xi − xj|| is the
distance between the points xi and xj, and w(xi, xj) is an edge correction weight. Here, we used the translation edge
correction.36
The empty space function F(r) is related to the probability that an arbitrary location s in the observation window has
an empty disk of radius r, b(s, r), around it. It is defined as
F(r) = 1 − P(the number of xi in b(s, r) = 0).
The empty space function can be estimated using a number of sample points in the observation windowW . Typically, a
grid is used. Letm be the number of sample points and di the distance from the sample point i to its nearest point in the
point pattern. Then, an estimator for the empty space function is
F̂(r) = the number of di < rm .
This estimator is hampered with edge effects since we cannot observe if a disk close to the boundary of the observation
windowwould have any points of the process outside thewindow. Therefore, we used theKaplan-Meiermethod to correct
for the edge effects.37
2.4 Descriptive statistics of the point patterns
We first estimated the pair correlation function for each of the sweat gland patterns shown in Figure 3 and thereafter,
obtained the groupwise pair correlation functions (see Figure 4) by pooling the estimated pair correlation functions of
all the subjects within the group by taking a weighted average of the values of the pair correlation fiunctions at each
distance.35 Since the number of points in the point patterns within each group varied quite a lot, the individual pair
correlation functions were weighted by the squared number of points, not number of points, when pooling.38-40 The pair
KURONEN et al. 7
F I GURE 4 Pooled pair correlation functions for the three
groups [Color figure can be viewed at wileyonlinelibrary.com]
0.00
0.25
0.50
0.75
1.00
0 100 200 300 400 500
r
g
(r
)
Controls MNA MNA Diagnosed
correlation functions show a clear sign of inhibition in all three groups (g(r)< 1 for small r). Further, the initial peak of
the functions appears approximately at the same distance for the control and suspected neuropathy groups. However,
the diagnosed neuropathy group has the initial peak at a slightly longer distance, indicating somewhat larger range of
inhibition than in the other two groups.
At very short distances, especially the control subjects seem to have some unexpected close pairs of points. Upon a
closer inspection of the point patterns and the videos, it was reasonable to assume that some sweat spots had been detected
as two nearby spots, instead of havingmerged into one. An obvious, simple idea to remove such close pairs of spots would
be tomerge all small glands having a larger spot closer than at some specified distancewith the larger spot. However, such
erroneous pairs of glands may appear at various (small) distances from each other and thus, applying a global limiting
distance is not reasonable. Instead of using an additional image analysis step, we include some of this inaccuracy in the
modeling and/or parameter estimation.
3 SEQUENTIAL POINT PROCESS MODEL
Since sweat glands activate at different times, we modeled the activation by using sequential point processes similar to
those suggested for modeling eyemovements.28 The points, activated sweat glands in our case, are generated sequentially
conditioning on the already existing points. Points are added until the observed number of points in the pattern has been
reached and the main focus here is to make inference on the arrival density. Below, we first recall the general sequential
model28 (Section 3.1) and specify it in our case without (Section 3.2) and with noise (3.4). Further, we discuss efficient
inference for the sequential models (Section 3.3) and, finally, fit the sequential model with noise to the sweat gland data
(Section 3.5).
3.1 General sequential point process model
Denote byW the observation window and by n the fixed number of points in the pattern. The first point x1 is assumed
to be uniformly distributed in W and the kth point, k= 2, … ,n, is assumed to follow the density y → f (y; x⃗k−1), where
x⃗k−1 = (x1, x2, … , xk−1). The density function for the whole point pattern (x1, … , xN) is then
x⃗n → 1|W |
n∏
k=2
f (xk; x⃗k−1),
8 KURONEN et al.
where 1/|W | is the contribution of the first point. A nice feature of the sequential point process models is that they have
a tractable likelihood even though it can be costly to compute.28
3.2 Soft-core model
The function f above should be chosen based on the phenomenon we would like to model. The activated sweat gland
location patterns are repulsive. Our first attempt was to use a hard-core model, where sweat glands cannot be closer
together than some minimum hard-core distance, but such a model turned out not to be flexible enough. Therefore, we
suggest to use a soft-core model with the density
fSC(y; x⃗k,R, 𝜅) ∝ exp
(
−
k∑
i=1
(
R
d(y, xi)
)2∕𝜅)
for adding the point y in the realization. Above, R> 0 is an inhibition range parameter and 0 < 𝜅 < 1 in the exponent
describes how “soft-core” the model is. In the limit as 𝜅 → 0, we obtain a hard-core process with hard-core distance
R. Some soft-core Gibbs point process models have been introduced in the literature,41,42 including models with the
particular interaction function that we use here.43
The log likelihood of the model becomes
l(R, 𝜅; x⃗n) = −log |W | − n∑
k=2
k−1∑
i=1
(
R
d(xk, xi)
)2∕𝜅
−
n∑
k=2
logZ(R, 𝜅, x⃗k), (1)
where
Z(R, 𝜅, x⃗k)−1 = ∫W exp
(
−
k−1∑
i=1
(
R
d(y, xi)
)2∕𝜅)
dy
is a normalizing constant.
3.3 Efficient likelihood inference for the sequential models
Even though the likelihood of a sequential point process can be costly to compute, the particular sum structure
in (1) allows faster computations. Using an integration scheme with J integration points y1, y2, … , yJ with weights
w1,w2, … ,wJ , the last term in (1) can be written as
n∑
k=2
logZ(R, 𝜅, x⃗k) =
n∑
k=2
log
(
∫W exp
(
−
k−1∑
i=1
(
R
d(y, xi)
)2∕𝜅)
dy
)−1
= −
n∑
k=2
log
J∑
j=1
wj exp
(
−
k−1∑
i=1
(
R
d(yj, xi)
)2∕𝜅)
.
In total, there are Jn(n− 1)/2 summands, among which only Jn are distinct. Therefore, the integrals are efficiently
calculated by evaluating the terms in the innermost sum only once.
3.4 Soft-core model with noise
To account for the incorrectly identified close pairs in the extracted point patterns, we used amixture model where one of
the components is a uniformly distributed error component. Such an error component can be added to any point process
KURONEN et al. 9
F I GURE 5 Parameter
estimates of the soft-core
model without and with noise
fitted separately to each subject
of the three groups (subject
numbers shown on the left)
R
Controls
MNA
MNA Diagnosed
20 40 60 0.2 0.4 0.6 0.8 0.0 0.1 0.2 0.3 0.4
205
203
149
97
96
71
61
40
20
10
73
50
42
36
23
With noise Without noise
model and here, we add it in the sequential soft-coremodel. The arrival density of a point y (after the uniformly distributed
first point) is then
fM(y; x⃗k,R, 𝜅, 𝜃) = (1 − 𝜃)fSC(y; x⃗k,R, 𝜅) + 𝜃|W |
= (1 − 𝜃)Z(R, 𝜅, x⃗k) exp
(
−
k∑
i=1
(
R
d(y, xi)
)2∕𝜅)
+ 𝜃|W | .
Therefore, the point at y comes from the soft-core process with probability 1 − 𝜃 (the first term on the right-hand side
of the formula) and from the uniformly distributed error process with probability 𝜃. Even though this model allows extra
points everywhere, not only near the real activated glands, it can improve estimation of the parameters as shown below.
However, the parameter 𝜃 cannot be interpreted directly as the probability of incorrectly identified glands since some of
the points without close neighbors regarded as noise could as well be true glands.
The log-likelihood of the soft-core model with uniformly distributed error is given by
lM(R, 𝜅, 𝜃; x⃗n) = − log |W | + n∑
k=2
log fM(xk; x⃗k−1,R, 𝜅, 𝜃). (2)
3.5 Application to the sweat gland data
The soft-coremodelwas fittedwithout andwithnoise to each sweat glandpoint pattern independently. First, we compared
the maximum likelihood estimates of the soft-core parameters obtained without or with added noise. Then, we fitted the
model with noise to the data in a Bayesian framework to be able to better compare the goodness-of fit of the sequential
soft-core model and the generative model presented in Section 4. We used regular grid based integration with 10 800
integration points to evaluate the likelihood in all cases.
3.5.1 Parameter estimates without and with added noise
The parameter estimates obtained by maximizing the log likelihood (1) or (2) with respect to the parameters can be seen
in Figure 5, where circles belong to the sequential soft-core model without noise and the pluses to the model with noise.
The estimates obtained without noise for the range parameter R are on average smaller and the “softness” parameter 𝜅
larger in the control group than in the neuropathy groups. However, for the model fitted with noise (circles on Figure 5),
only the mixture parameter 𝜃, which is estimated larger for the control group than for the neuropathy groups, differs
between the groups. Note that the 𝜅 and R parameters coincide for the twomodels if the mixture parameter 𝜃 is estimated
to be zero, see subjects 42, 50, and 71 in Figure 5.
10 KURONEN et al.
Without noise With noise
0 100 200 300 0 100 200 300
0.0
0.5
1.0
1.5
r
g
(r
)
F IGURE 6 Empirical pair
correlation functions (black
lines) for subject 205 in the end
of the video recording together
with 95% global envelopes (grey
areas) constructed from 25 000
simulations from the soft-core
model estimated without (left)
and with (right) noise
We investigated the goodness-of-fit of the fitted soft-core models by using the pair-correlation function. We generated
samples from the sequential soft-core models with parameters R and 𝜅 estimated with and without noise. The uniform
noise was not simulated. Figure 6 shows the empirical pair-correlation functions for subject 205 for the soft-core model
estimatedwith andwithout noise together with 95% global envelopes44,45 calculated from 25 000 samples of eachmodel. It
can be seen that for this subject, the range parameter is clearly underestimated if estimation is donewithout accounting for
noise. For the other subjects, the goodness-of-fit of themodel with noise was also as good or better than the goodness-of-fit
of the model without noise. The difference between the observed and model (fitted with noise) based pair correlation
function at short distances is explained by the incorrectly recorded close pairs of points that are present in the data but
not in the simulations.
3.5.2 Bayesian inference of the model with noise
We fitted the soft-coremodelwith noise to the sweat gland data also by using standard likelihood-basedBayesian approach
with robust adaptive Metropolis algorithm.46 We ran the MCMC for 120 000 iterations and discarded the first 20 000
iterations as burn-in. As the prior distribution for the range parameter R, we used the Gamma distribution with shape
parameter 3 and scale parameter 70/3 and the priors for 𝜅 and 𝜃were both the uniform distribution on [0, 1]. The posterior
histograms in Figure 7 show some variation within the groups but no clear differences between the groups: The arrival
density parameters R and 𝜅 were estimated to be rather similar in all groups. The 𝜃 parameter related to the errors appears
to be somewhat larger in the control group than in the other two groups.
Figure 8 shows the empirical pair-correlation functions for each subject together with the global envelopes44,47 cal-
culated from 25 000 simulations from the posterior predictive distribution of the fitted soft-core models with noise.
In most cases, the envelopes cover the empirical curves. For some subjects, especially for the controls, the empirical
pair-correlation function is not covered by the envelopes at very short distances. This is expected, as mentioned earlier,
since according to the model used, this behavior is caused mainly by noise, which was not simulated. Thus, the pair cor-
relation function reveals that, ignoring the noise at short distances, the quite clear inhibition and the gradual increase of
interpoint distances around 50-75 pixels in the empirical patterns is captured by the soft-core model. However, the char-
acteristic peak in the pair correlation functions around 100 pixels is not always captured by the model, particularly for
the patients who have smaller number of activated sweat glands, as the envelopes are quite wide close to the peak. All the
empirical pair correlation functions have the same characteristic behavior.
4 GENERATIVE POINT PROCESS MODEL
In our second approach, we first model the underlying unobserved sweat glands and then, model the activated sweat
glands as an independent thinning of the underlying gland pattern. Modeling the glands and the activation of them
separately allows one to answer questions regarding specifically the activation process. One possible hypothesis is that the
KURONEN et al. 11
R
96 Controls
97 Controls
149 Controls
203 Controls
205 Controls
10 MNA
20 MNA
40 MNA
61 MNA
71 MNA
23 MNA Diagnosed
36 MNA Diagnosed
42 MNA Diagnosed
50 MNA Diagnosed
73 MNA Diagnosed
Prior
25 50 75 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00
F IGURE 7 Posterior marginals for each subject (subject number and group given on the right) and each parameter (R, 𝜅, and 𝜃) for the
soft-core model estimated with noise. The prior densities are given on the last row
MNA Diagnosed, 23 MNA Diagnosed, 36 MNA Diagnosed, 42 MNA Diagnosed, 50 MNA Diagnosed, 73
MNA, 10 MNA, 20 MNA, 40 MNA, 61 MNA, 71
Controls, 96 Controls, 97 Controls, 149 Controls, 203 Controls, 205
0 100 200 300 400 0 100 200 300 400 0 100 200 300 400 0 100 200 300 400 0 100 200 300 400
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
r [px]
g
(r
)
F IGURE 8 Empirical pair correlation functions (black lines) for each subject (subject number shown in the label above each plot) in
the end of the video recording together with 95% global envelopes (gray areas) constructed from 25 000 simulations from the posterior
predictive distribution of the soft-core model estimated with noise
underlying gland pattern itself is not different between controls and subjects with neuropathy, but the activation process
is different. More specifically, almost all glands should activate on healthy subjects while the glands of the subjects with
neuropathy could have a tendency to leave larger holes in the activation process.16
4.1 Model specification
It seems reasonable to assume that the underlying (unobserved) sweat gland pattern is a rather densely packed reg-
ular point pattern covering the whole skin. To obtain such a structure, some type of soft-core sequential inhibition
process, where points are added as long as it is possible (we do not know the actual number of glands), would be
12 KURONEN et al.
appropriate. However, it is not straightforward to decide when to stop adding points since theoretically, soft-core type
of interaction always allows new points. Instead, we start by generating a simple sequential inhibition (SSI) model,35
which is then disturbed to obtain a soft-core structure. A sample from the SSI model is generated sequentially by
proposing points from the uniform distribution and accepting them if the pattern satisfies the hard-core condition with
hard-core distance R, that is, the new proposed point does not lie within distance R from any earlier point. This is
continued until there is no space left for new points. The disturbed SSI model is obtained from the “pure” SSI model
by displacing the location of each point with an independent zero mean isotropic Gaussian random variable with
covariance 𝜎2I.
We assume independent gland activation, that is, that the final pattern is a result of an independent thinning of the
underlying disturbed SSI process. Therefore, the model has three parameters: inhibition range R, hardness of inhibition
𝜎, and probability of activation p.
4.2 Parameter estimation using approximate Bayesian computation
For the generative model, we cannot write down the likelihood. However, sampling from the model is easy. We used
the method proposed by Wang48 to generate samples from the SSI process. Approximate Bayesian computation (ABC)
is a method for Bayesian inference in situations where the likelihood of the model is intractable,49,50 but it is possible to
simulate the model. It is based on sampling from the (pseudo-) posterior distribution
𝜋𝜖(𝜃) = 𝜋(𝜃)P(||s(Y𝜃) − s(y)|| < 𝜖),
where Y𝜃 follows the model with parameter vector 𝜃, y is the data, 𝜋(⋅) is the prior distribution for the parameters, s is an
appropriately chosen summary statistic, and 𝜖 is a tolerance level.
4.2.1 ABC-MCMC
AsimpleABC rejection sampler is expressed inAlgorithm1. This basic algorithmcan be rather inefficient, but fortunately,
there are severalmore efficient algorithms for performingABC.Weused an adaptiveABC-MCMCalgorithm.51 In our data
study below, the MCMC was run for 10 000 000 iterations and the 250 000 simulated parameter values with the smallest
distances ||s(z)− s(y)|| were taken as the posterior sample.
Algorithm 1. A simple ABC rejection sampler
for i ← 1,M do
repeat
Generate parameter vector 𝜃′ from the prior distribution 𝜋
Generate a realization z from the model with parameter vector 𝜃′
until ‖s(z) − s(y)‖ ≤ 𝜖
𝜃i ← 𝜃′
end for
4.2.2 Summary statistics
The choice of summary statistics is crucial for the ABC method to work. For a regular point process model, it is natural
to use summary statistics based on the pair correlation function g. Instead of using the full pair correlation function g, we
tried to find a specific part of it that would be sufficient for our purpose following the rule of thumb52 that the number of
summary statistics in theABC approach should approximatelymatchwith the number of parameters to be estimated. The
location of the first peak of the pair correlation function is intuitively connected to the inhibition range R. However, the
KURONEN et al. 13
R p
96 Controls
97 Controls
149 Controls
203 Controls
205 Controls
10 MNA
20 MNA
40 MNA
61 MNA
71 MNA
23 MNA Diagnosed
36 MNA Diagnosed
42 MNA Diagnosed
50 MNA Diagnosed
73 MNA Diagnosed
Prior
40 60 80 100 0.00 0.25 0.50 0.75 1.00 0 10 20 30 40
F IGURE 9 Posterior marginals for each subject (subject number and group given on the right) and each parameter (R, p, and 𝜎) of the
generative model
location of the first peak can be difficult to estimate exactly and thus, we used the smallest distance r1 > 10 pixels where
g(r1)= 0.75 as the location of the uphill before the first peak. Furthermore, the slope of the uphill provides information
on the “softness” parameter 𝜎 and we chose the smallest distance r2 > 10 pixels where g(r2)= 1 as the second summary
statistic. Finally, the smallest distance r3 in the empty space function F where F(r3)= 0.5 was taken as the third summary
statistic to represent the activation probability p. The empty space functionwas chosen because it gives information on the
number of points but is not greatly affected by erroneous nearby points. Since all the chosen summary statistics, r1, r2, and
r3, have a similar order of magnitude, we did not have to add any weights in the ABC algorithm. The specific values 0.75
and 1 were chosen to be somewhat separated and not too small to account for possible errors caused by splitting of spots
into multiple glands that would cause the pair-correlation function not to start from zero. In addition, we only considered
distances greater than 10 pixels since at very short distances the kernel estimator of the pair correlation function is not
very reliable. These choices worked well for the sweat gland data, as demonstrated below.
4.3 Application to the sweat gland data
The generativemodel was fitted to the sweat gland data using the ABC approach described above. In addition to the above
specifications, we needed to set the priors. For R we used an improper, uniform prior on [40,∞) restricting that R could
not be unreasonably small, while in addition to being unrealistic, small R values result in a large number of points in the
SSI process which is computationally challenging. The prior of p was uniform on [0.1, 1], stating that at least 10% of the
glands (modeled by the underlying disturbed SSI process) needed to activate and thus be observed. Furthermore, for 𝜎,
we used the gamma distribution with the shape parameter equal to 10/3 and scale parameter equal to 3. While the priors
R and p can be considered rather non-informative, the prior for 𝜎 was somewhat informative suggesting positive, but not
too large 𝜎. Note that if 𝜎 was very large in comparison to R, it would break all the structure of the SSI process, which is
unreasonable.
The posterior marginal histograms for the parameters can be seen in Figure 9 and 95% global envelopes for the pair
correlation function constructed from 25 000 simulations from the posterior predictive distribution in Figure 10. As can
be seen in Figure 9, the parameter estimates vary somewhat between the subjects and groups. Differences in the softness
of the model, that is, in the values of the parameter 𝜎, are small. However, there seems to be a slight tendency for the
inhibition range R to be a little smaller in the control group than in the MNA groups, but the difference is not clear based
on the limited amount of data we have. The range was always between 60 pixels and 100 pixels. Furthermore, the control
subjects tend to have a larger activation probability than theMNA patients, but the within group variation is large. This is
14 KURONEN et al.
MNA Diagnosed, 23 MNA Diagnosed, 36 MNA Diagnosed, 42 MNA Diagnosed, 50 MNA Diagnosed, 73
MNA, 10 MNA, 20 MNA, 40 MNA, 61 MNA, 71
Controls, 96 Controls, 97 Controls, 149 Controls, 203 Controls, 205
0 100 200 300 400 0 100 200 300 400 0 100 200 300 400 0 100 200 300 400 0 100 200 300 400
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
0.0
0.5
1.0
1.5
r [px]
g
(r
)
F IGURE 10 Empirical pair correlation functions (lines) for each subject (subject number shown in the label above each plot) in the
end of the video recording together with 95% global envelopes (gray areas) and means (dashed lines) constructed from 25 000 simulations
from the posterior predictive distribution of the generative model
in agreement with earlier studies, which indicate that a larger number of sweat glands of controls than of MNA patients
activate.4
According to the visual evaluation of the global envelopes of the pair-correlation function (see Figure 10) and empty
space function (see Figure 11), the model seems to fit quite well to the data. It captures the behavior of the pair cor-
relation function both at small distances and around the initial peak. It should also be mentioned that the envelopes
for the pair correlation function are rather wide at small distances covering the observed functions almost in all cases,
even though the model did not include any error term. The wide envelopes are due to the relatively wide posterior
distribution of 𝜎. Namely, large 𝜎 can lead to some close pairs in the patterns and consequently also positive values
of the pair correlation function at small distances. Another reason for the relatively wide envelopes may be that the
summary statistics used in the ABC approach were chosen such that they do not use any information at very short
distances.
We explored a few other priors for 𝜎, namely, improper uniform and exponential distributions with means 1, 2, and 4.
The posterior distributions of the other parameters were not affected by the choice of the prior for 𝜎, but the posterior of 𝜎
itself was somewhat sensitive to the choice and also the goodness-of-fit of themodel measured by the pair correlation was
affected. Namely, the improper uniform prior led to wider posterior distribution of 𝜎 and large 𝜎 caused the variation of
the pair correlation function to be even higher at small distances. On the other hand, the strict exponential priors shrank
the posterior distribution toward zero, and very small 𝜎 caused the peak of the pair correlation function to be too sharp.
Thus, the disturbance parameter 𝜎 needed a somewhat informative prior to lead to a good fit of the model.
We simulated patterns from the posterior predictive distribution and the simulated patterns mimic the data patterns
rather well, see Figure 12. Note, in particular, that the independent thinning seems to produce empty regions similar to
those observed in the point pattern data, as also indicated by the empty space function (Figure 11).
5 DISCUSSION
Videos produced by dynamic sweat tests provide not only the total and per gland sweat volume and the number and
density of activated sweat glands but also a time series of spatial patterns of the activated sweat glands. Up to now, only
KURONEN et al. 15
MNA Diagnosed, 23 MNA Diagnosed, 36 MNA Diagnosed, 42 MNA Diagnosed, 50 MNA Diagnosed, 73
MNA, 10 MNA, 20 MNA, 40 MNA, 61 MNA, 71
Controls, 96 Controls, 97 Controls, 149 Controls, 203 Controls, 205
0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
r [px]
F
(r
)
F IGURE 11 Empirical empty space functions (black lines) for each subject (subject number shown in label) in the end of the video
recording together with 95% global envelopes (gray areas) and means (dashed lines) constructed from 25 000 simulations from the posterior
predictive distribution of the generative model
the sweat volume and the number/density of activated sweat glands have been used in the analysis and comparison
between healthy subjects and subjects with neuropathy. However, visual inspection has indicated that the spatial pat-
tern of activated sweat glands of subjects with neuropathy seem to have more empty areas than the patterns of healthy
controls. This may indicate that the neuropathy does not result in random deactivation of sweat glands but the glands
that activate are arranged in a different pattern than in control patterns. To quantify this observation, spatial analysis is
needed.
Here, we analyzed videos of sweat gland activation recorded from five controls, five subjects with suspected neuropa-
thy, and five subjects with diagnosed neuropathy. The initial non-parametric spatial analysis by using pair correlation
functions showed clear indication of inhibition in all three groups. In addition, it suggested some differences between the
spatial patterns of activated sweat glands in subjects diagnosed for neuropathy and the non-diagnosed and control groups.
To further investigate and compare the spatial patterns and the activation processes, we suggested two point process
models for the activation of sweat glands, a sequential soft-core model describing the appearance of the activated sweat
glands and a thinned disturbed SSI process, that we call a generative model, where we start by modeling the underlying
unobserved sweat gland pattern.
Maximizing the log-likelihood function of a sequential point process has been considered computationally costly due
to the integrals in the normalizing constants.28 However, for the sequential soft-coremodel, these integrals have a particu-
lar sum form which allows efficient computation of the log-likelihood and in turn Bayesian inference. The same efficient
computation scheme is applicable for any sequential point process having an arrival density with a similar sum structure.
To estimate the parameters of the generative model, we employed an ABC algorithm since the likelihood function was
not easily available.
Even though our proposed image analysis approach worked well, there were some incorrectly identified close pairs of
glands in the extracted point patterns. To take into account such errors, we added an error term in the sequential soft-core
model resulting in a mixture model having a soft-core component and a uniform noise component. For the generative
model, on the other hand, the summary statistics in the ABC approach were chosen such that they were robust to close
pairs of points. The good fit of the models shows that it is possible to account for some, but likely not many, errors in the
pre-processing of videos. Such an analysis needs to be applied with care, but allows the researcher to focus on the analysis
more than on the pre-processing.
16 KURONEN et al.
Controls, 96 MNA, 20 MNA Diagnosed, 42
data
m
odel
0 1000 2000 0 1000 2000 0 1000 2000
0
500
1000
1500
2000
0
500
1000
1500
2000
x [px]
y
[p
x]
F IGURE 12 The original point patterns (top) and patterns generated from the corresponding posterior predictive distributions of the
generative model (bottom) for one subject from each group (96, 20, and 42)
The proposed models were fitted to the data. The fit of the generative model was good regardless of the disease status,
as the point patterns in Figure 12 and the spatial summary characteristics in Figures 10 and 11 show. It should be noted
that this model based on independent activation of underlying sweat glands was also able to produce similar holes in the
point patterns as observed earlier in the sweat gland patterns of diabetic patients;16 see, for example, the bottom right
plot in Figure 12. Apparently the SSI process with random displacements of points is a working model for the underlying
sweat glands.
Themodel parameters estimated from the patterns of healthy subjects and of subjects suffering from neuropathy were
compared. Based on the generative model, detailed information was obtained on the sweat patterns: while the hardness
of inhibition parameter was rather similar for each subject, the inhibition range R and the activation probability p were
characteristic for each subject. The former characterizes the pattern of all glands, while the latter is the probability of
these glands to activate. There was some indication that the inhibition range R was slightly smaller and the activation
probability p larger in the control group than in the neuropathy groups (Figure 9). Both of these observations are in
agreement with an earlier study,16 where the density of activated sweat glands was found to be lower for diabetic patients
than controls. However, the spatial analysis provides more detailed information than the non-spatial analysis alone, and
the model parameters might facilitate classification of subjects in the future. For example, the control subject 203 had
rather small number of activated sweat glands (Figure 3), but simultaneously a small inhibition range parameter. Thus,
the parameters may contain more valuable information jointly as either of them alone.
We believe that themodels suggested here, especially the generativemodel, would givemore impactful insight into the
sweat patterns in a larger study includingmore subjects. Especially, combining the point pattern approach presented here
with non-spatial covariates like age or body mass index would provide a better understanding of the sweat patterns and
provide an explanation for the individually varying generative model parameters. In addition, incorporating the amount
of sweat produced over time by individual sweat glands into the analysis would further improve our understanding and
lead to further methodological development.
ACKNOWLEDGEMENTS
Mikko Kuronen and Mari Myllymäki have been financially supported by the Academy of Finland (Project Numbers
295100, 306875, and 327211) and Aila Särkkä by the Swedish Research Council (VR 2013-5212) and by the Swedish
KURONEN et al. 17
Foundation for Strategic Research (SSF AM13-0066). The authors thank Matti Vihola for useful discussions. We are also
grateful for the two anonymous reviewers for the valuable comments.
DATA AVAILABILITY STATEMENT
The original video data are not publicly available. The point pattern data and code used in the analysis is available as a
github repository: https://github.com/mikkoku/SweatPaper.git.
ORCID
Mikko Kuronen https://orcid.org/0000-0002-8089-7895
REFERENCES
1. Coon JM, Rothman S. The sweat response to drugs with nicotine-like action. J Pharmacol Exp Ther. 1941;73(1):1-11.
2. Lader M, Montagu J. The psycho-galvanic reflex: a pharmacological study of the peripheral mechanism. J Neurol Neurosurg Psychiatry.
1962;25:126-133.
3. Minor V. Ein neues Verfahren zu der klinischen Untersuchung der Schweißabsonderung. Deutsche Zeitschrift für Nervenheilkunde.
1928;101:302-307.
4. Loavenbruck AJ, Hodges JS, Provitera V, Nolano M, Wendelshafer-Crabb G, Kennedy WR. A device to measure secretion of individual
sweat glands for diagnosis of peripheral neuropathy. J Peripher Nerv Syst. 2017;22(2):139-148. https://doi.org/10.1111/jns.12212.
5. Sumner C, Sheth S, Griffin J, Cornblath D, Polydefkis M. The spectrum of neuropathy in diabetes and impaired glucose tolerance.
Neurology. 2003;60(1):108-111. https://doi.org/10.1212/WNL.60.1.108.
6. Said G. Diabetic neuropathy–a review. Nat Clin Pract Neurol. 2007;3:331-340. https://doi.org/10.1038/ncpneuro0504.
7. Sharma S, Venkitaraman R, Vas PRJ, Rayman G. Assessment of chemotherapy-induced peripheral neuropathy using the
LDIFLARE technique: a novel technique to detect neural small fiber dysfunction. Brain Behav. 2015;5(7):e00354. https://doi.org/10.1002/
brb3.354.
8. Hilz MJ, Dütsch M. Quantitative studies of autonomic function. Muscle Nerve Offic J Am Assoc Electrodiagnostic Med. 2006;33(1):6-20.
https://doi.org/10.1002/mus.20365.
9. Illigens BMW, Gibbons CH. Sweat testing to evaluate autonomic function. Clin Auton Res. 2009;19:79-87. https://doi.org/10.1007/s10286-
008-0506-8.
10. Minota K, Coon EA, Benarroch EE. Neurologic aspects of sweating and its disorders.Neurology. 2019;92(21):999-1005. https://doi.org/10.
1212/WNL.0000000000007540.
11. Fealey RD, Low PA, Thomas JE. Thermoregulatory sweating abnormalities in diabetes mellitus. Mayo Clin Proc. 1989;64(6):617-628.
https://doi.org/10.1016/S0025-6196(12)65338-5.
12. Papanas N, Papatheodorou K, Christakidis D, et al. Evaluation of a new indicator test for sudomotor function (Neuropad®) in the diagno-
sis of peripheral neuropathy in type 2 diabetic patients. Exp Clin Endocrinol Diabetes. 2005;113:195-198. https://doi.org/10.1055/s-2005-
837735.
13. Harris DR, Polk BF, Willis I. Evaluating sweat gland activity with imprint techniques. J Investig Dermatol. 1972;58(2):78-84. https://doi.
org/10.1111/1523-1747.ep12551676.
14. Low PA, Caskey PE, Tuck RR, Fealey RD, Dyck PJ. Quantitative sudomotor axon reflex test in normal and neuropathic subjects. Ann
Neurol. 1983;14:573-580.
15. Low VA, Sandroni P, Fealey RD, Low PA. Detection of small-fiber neuropathy by sudomotor testing. Muscle Nerve. 2006;34(1):57-61.
https://doi.org/10.1002/mus.20551.
16. Provitera V, Nolano M, Caporaso G, Stancanelli A, Santoro L, Kennedy WR. Evaluation of sudomotor function in diabetes using the
dynamic sweat test. Neurology. 2010;74(1):50-56. https://doi.org/10.1212/WNL.0b013e3181c7da4b.
17. KennedyWR, SelimMM,Wendelschafer-Crabb G, et al. A device to quantify sweat in single sweat glands to diagnose neuropathy. Journal
of Medical Devices. 2013;7(3):030941. https://doi.org/10.1115/1.4024527.
18. Loavenbruck A, Sit N, Provitera V, Kennedy W. High-resolution axon reflex sweat testing for diagnosis of neuropathy. Clin Auton Res.
2019;29:55-62. https://doi.org/10.1007/s10286-018-0546-7.
19. Byers S, Raftery AE. Nearest-neighbor clutter removal for estimating features in spatial point processes. J Am Stat Assoc.
1998;93(442):577-584. https://doi.org/10.1080/01621459.1998.10473711.
20. Cressie N, Lawson AB. Hierarchical probability models and Bayesian analysis of mine locations. Adv Appl Probab. 2000;32(02):315-330.
https://doi.org/10.1017/s0001867800009940.
21. Walsh DCI, Raftery AE. Detecting mines in minefields with linear characteristics. Technometrics. 2002;44(1):34-44.
22. Redenbach C, Särkkä A, Sormani M. Classification of points in superpositions of Strauss and Poisson processes. Spat Stat. 2015;12:81-95.
https://doi.org/10.1016/j.spasta.2015.03.003.
23. Rajala T, Redenbach C, Särkkä A, Sormani M. Variational Bayes approach for classification of points in superpositions of point processes.
Spat Stat. 2016;15:85-99. https://doi.org/10.1016/j.spasta.2015.12.001.
24. Lund J, Rudemo M. Models for point processes observed with noise. Biometrika. 2000;87:235-249. https://doi.org/10.1093/biomet/87.2.
235.
18 KURONEN et al.
25. Lund J, Penttinen A, Rudemo M. Bayesian analysis of spatial point patterns from noisy observations. [Ph.D. thesis]. Statistical inference
and perfect simulation for point processes observed with noise by J Lund. Denmark: Department of Mathematics and Physics, Royal
Veterinary and Agricultural University; 1999.
26. Chakraborty A, Gelfand AE. Analyzing spatial point patterns subject to measurement error. Bayesian Anal. 2010;5(1):97-122. https://doi.
org/10.1214/10-BA504.
27. Bar-Hen A, Chadœuf J, Dessard H, Monestiez P. Estimating second order characteristics of point processes with known independent
noise. Stat Comput. 2013;23(3):297-309. https://doi.org/10.1007/s11222-011-9311-7.
28. Penttinen A, Ylitalo AK. Deducing self-interaction in eye movement data using sequential spatial point processes. Spat Stat. 2016;17:1-21.
https://doi.org/10.1016/j.spasta.2016.03.005.
29. Bezanson J, Edelman A, Karpinski S, Shah VB. Julia: a fresh approach to numerical computing. SIAM Rev. 2017;59(1):65-98. https://doi.
org/10.1137/141000671.
30. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2018.
31. Hinkley DV. Inference about the change-point in a sequence of random variables. Biometrika. 1970;57(1):1-17. https://doi.org/10.2307/
2334932.
32. Yao YC, Davis RA. The asymptotic behavior of the likelihood ratio statistic for testing a shift in mean in a sequence of independent normal
variates. Sankhya Ind J Stat Ser A (1961-2002). 1986;48(3):339-353.
33. James B, James KL, Siegmund D. Tests for a change-point. Biometrika. 1987;74(1):71-83. https://doi.org/10.2307/2336022.
34. Pettitt AN. A non-parametric approach to the change-point problem. J Royal Stat Soc Ser C (Appl Stat). 1979;28(2):126-135. https://doi.
org/10.2307/2346729.
35. Illian J, Penttinen A, Stoyan H, Stoyan D. Statistical Analysis and Modelling of Spatial Point Patterns. Chichester, England: John Wiley &
Sons, Ltd; 2008.
36. Stoyan D, Fractals SH. Random Shapes and Point Fields: Methods of Geometrical Statistics. Hoboken, NJ: Wiley; 1994.
37. Baddeley A, Gill RD. Kaplan-Meier estimators of distance distributions for spatial point processes. Ann Stat. 1997;25(1):263-292.
38. Diggle PJ, Mateu J, Clough HE. A comparison between parametric and non-parametric approaches to the analysis of replicated spatial
point patterns. Adv Appl Probab. 2000;32:331-343.
39. Myllymäki M, Panoutsopoulou IG, Särkkä A. Analysis of spatial structure of epidermal nerve entry point patterns based on replicated
data. J Microsc. 2012;247(3):228-239.
40. Schladitz K, Särkkä A, Pavenstädt I, Haferkamp O, Mattfeldt T. Statistical analysis of intramembranous particles using freeze fracture
specimens. J Microsc. 2003;211:137-153.
41. Ogata Y, TanemuraM. Estimation of interaction potentials of spatial point patterns through the maximum likelihood procedure. Ann Inst
Stat Math. 1981;33(2):315-338. https://doi.org/10.1007/BF02480944.
42. Ogata Y, Tanemura M. Likelihood analysis of spatial point patterns. J Royal Stat Soc Ser B (Methodol). 1984;46(3):496-518. https://doi.org/
10.1111/j.2517-6161.1984.tb01322.x.
43. Baddeley A, Rubak E, Turner R. Spatial Point Patterns: Methodology and Applications with R. London, UK: Chapman & Hall/CRC Press;
2015.
44. MyllymäkiM,Mrkvicˇka T, Seijo H, Grabarnik P, HahnU. Global envelope tests for spatial processes. J Royal Stat Soc Ser B (StatMethodol).
2017;79:381-404. https://doi.org/10.1111/rssb.12172.
45. Myllymäki M, Mrkvicˇka T. GET: global envelopes in R; 2019. arXiv:1911.06583 [stat.ME].
46. Vihola M. Robust adaptive Metropolis algorithm with coerced acceptance rate. Stat Comput. 2012;22:997-1008. https://doi.org/10.1007/
s11222-011-9269-5.
47. Mrkvicˇka T, Myllymäki M, Jílek M, Hahn U. A one-way ANOVA test for functional data with graphical interpretation; 2018.
arXiv:1506.01646 [stat.ME].
48. Wang JS. A fast algorithm for random sequential adsorption of discs. Int J Modern Phys C. 1994;05(04):707-715. https://doi.org/10.1142/
S0129183194000817.
49. Marin JM, Pudlo P, Robert CP, Ryder RJ. Approximate Bayesian computational methods. Stat Comput. 2012;22(6):1167-1180. https://doi.
org/10.1007/s11222-011-9288-2.
50. Sunnåker M, Busetto AG, Numminen E, Corander J, Foll M, Dessimoz C. Approximate Bayesian computation. PLoS Comput Biol.
2013;9(1):1-10. https://doi.org/10.1371/journal.pcbi.1002803.
51. Vihola M, Franks J. On the use of approximate Bayesian computation Markov chain Monte Carlo with inflated tolerance and
post-correction. Biometrika. 2020;107(2):381-395. https://doi.org/10.1093/biomet/asz078.
52. Li W, Fearnhead P. On the asymptotic efficiency of approximate Bayesian computation estimators. Biometrika. 2018;105(2):285-299.
https://doi.org/10.1093/biomet/asx078.
How to cite this article: Kuronen M, Myllymäki M, Loavenbruck A, Särkkä A. Point process models for sweat
gland activation observed with noise. Statistics in Medicine. 2021;1–18. https://doi.org/10.1002/sim.8891