Introduction

A study was performed by Santello et al. in 1998 to investigate whether grasping postures showed natural patterns of coordination between the many joints of the human hand. Participants were asked to grasp imaginary objects of diverse shapes with one of their hands. The angles between various segments of the hand were measured with a tracking glove, and analyzed to identify ``postural synergies". This report documents a replication of the 1998 study, using visual tracking of the hand with a wearable headset. The key analysis of interest of the replication is a principal component analysis (PCA) of the joint angles data from grasping of a number of geometries. An important finding of the original paper was that the first two principal components account for more than 80% of the variance in posture.

While a goal of the original experiment was to better understand the neuromuscular basis of movement, its results have proved significant for efficient motion planning in robotic manipulators, as well as for the design and control of prosthetic hands. The choice of this study for replication was motivated by an interest in simplifying the control of prosthetic devices by exploiting an understanding of the neural control of movement.

It was anticipated that the use of visual tracking might lead to less accurate data for fingers which were partially or completely occluded from the camera. Since the stimuli were entirely verbal, it was likely that different subjects would interpret them differently based on their own imagined shape and size of each object. Finally, owing to a lack of experience with implementation, minor challenges were expected in obtaining joint angle data in a usable form from a proprietary headset, and performing a PCA in the chosen computing environment.

Project repository: https://github.com/psych251/santello1998

Original paper: https://github.com/psych251/santello1998/blob/main/original_paper/Santello1998_postural_synergies.pdf

Preregistration: https://osf.io/97yvz

Methods

Planned Sample

The sample of the original study as well as the replication was 5 right-handed participants. The age range of the original subjects was 30-41 years, while the replication used subjects of ages 22-30 years.

Materials

The original paper used a Cyberglove hand-tracking device while the replication used a Microsoft Hololens 2. A keyboard was used in place of a contact switch to register the completion of the grasping task.

Procedure

“Subjects were instructed to shape the right hand to grasp and use a large number of imagined objects (n = 57; Table 1). They were encouraged to move the proximal arm in tandem with the hand motion. An object was named at the beginning of each trial. We selected objects spanning a large range of sizes and shapes to assess the consequent modulation of hand posture…The subject was asked to imagine the object floating in space at a distance of ∼40 cm anterior to the subject’s frontal plane. The elbow and wrist rested on a flat surface, the forearm was horizontal, the arm was oriented in the parasagittal plane passing through the shoulder, and the hand was in a semipronated position. At presentation of a “go” signal, the subject moved the arm and hand as if to grasp and use the named object. A contact switch was released at movement onset. When subjects had attained a static hand posture, they pressed a second switch with their left hand. Each subject performed a total of five trials for each of the objects; all trials were presented in random order."

The above procedure was followed with the following exceptions:

  1. The subject was asked to wear a Microsoft Hololens 2 throughout the procedure, and ensure that their entire right hand was in the field of view before confirming task completion.
  2. No contact switch was used at movement onset.
  3. The spacebar on a keyboard was pressed with the left hand to indicate completion of the task.
  4. Three (not five) trials were presented for each object.

The list of objects presented in Table 1 of the paper was used without modification.

Analysis Plan

The key analysis of interest is the principal component analysis of the static joint angles after the goal pose is achieved. Hence, only one frame of joint angles data was used in the analysis for each trial.

The joint angles were recorded using a Unity app written for this purpose, with the hand-tracking feature built into the Windows Mixed Reality Toolkit (MRTK). A C# script was written to record joint positions when the participant presses the spacebar on a keyboard, compute joint angles from the positions, and log them to a text file as comma-separated values. The script is available at https://github.com/psych251/santello1998/blob/main/code/hand_tracking/Assets/logJointData.cs.

The data was placed into Microsoft Office Excel spreadsheets by the following procedure:

  1. Before data collection: populating the ParticipantID column manually, and randomizing a list of objects with three instances of each object in the second column (using the Excel RAND() function to assign a random numeric value to each entry, followed by sorting numbers and their accompanying objects in decreasing order); repeating each instance 15 times so that the joint angles can be pasted after data collection;
  2. After data collection: manually copying and pasting comma-separated data from the log file into the third column, using the Text Import Wizard to separate trial number, timestamp, joint name, and joint angle into separate columns; selecting and deleting alternate blank rows; replacing perfect zeros (representing failed data capture) with blank cells.

The joint angles estimated from visual tracking were (following the original work):

  1. metacarpal-phalangeal (mcp) and proximal interphalangeal (pip) flexion of each of the 4 fingers
  2. abduction between the 3 pairs of adjacent fingers
  3. mcp flexion, abduction, ip flexion, and rotation of the thumb

Abduction of the thumb was computed as the angle between (a) the metacarpal bone of the index finger, and (b) the projection of the metacarpal bone of the thumb into the plane of the palm. Rotation of the thumb was computed as the angle between the thumb metacarpus and the plane of the palm, defined as positive when the thumb is closer to the ventral side. The plane of the palm, for these two calculations, was defined by its normal vector, which in turn was defined to be perpendicular to both the index and middle finger metacarpal bones.

The following analysis procedure was followed precisely, except with three trials rather than five:

“For each subject, the five trials per object were first averaged to obtain a total of 57 hand postures. For each sensor, we then subtracted the grand mean computed over the 57 objects (so that the range was centered about 0°). The hand posture for each object was thus characterized as a “waveform” of the values of the 15 sensors (Santello and Soechting, 1997). The principal components (PCs) were then computed from the eigenvalues and eigenvectors of the matrix of the covariance coefficients between each of the 57 waveforms."

Differences from Original Study

The replication differed from the original study in that:

  1. data was collected using a headset rather than a tracking glove,
  2. participant ages belonged to a lower range than the original, and
  3. each object was presented three times and not five times.

The conclusions are expected to be robust to these changes.

Methods Addendum (Post Data Collection)

Verbal instructions to participants were as follows.

  1. Imagine the object is about 40 centimeters in front of you and you are using it with one hand.
  2. Try to keep the virtual information tab showing frame rate and CPU usage below your hand in your field-of-vision so that the chances of joint angles being captured is higher.
  3. If the information tab disappears at any time, stop and inform the experimenter.
  4. If you do not know what the object is, you can ask or imagine your closest guess.
  5. If you do not hear or understand what was said, please ask again.
  6. If you accidentally press the spacebar twice or miss pressing it, inform the experimenter so that it can be marked on the object list as a duplicate recording.

Participant 2 made a double-press twice, and participant 3 and 5 once each. The trial number was noted and matched the total number of readings in the log file. The experiment with participant 3 was interrupted by a draining of the Hololens battery and needed to be repeated from the beginning. The interruption recurred with participant 4; however, it was between two key presses, and therefore could be resumed without repetition.

The recordings corresponding to double-presses were excluded before analysis.

The original script for converting joint positions to joint angles calculated the three angles (pairwise) between the proximal phalanges of adjacent fingers. While these angles are valid measurements of hand posture and the result is expected to be robust to how they are measured, they are not ‘true’ abduction angles. It was thought that projecting the phalange vector onto the plane of the palm before measuring the angle might yield qualitative analyses closer to those in the original study where a glove was used to measure abduction angles. Hence, a measurement of the projected angle was added to the logging script in order to make exploratory analysis possible in this regard. This change was made before data collection from the fifth and final participant. It was made without any data processing or analysis of the previous four participants, and only additions were made; the sections of code used to make the pre-registered measurements were not edited in any way.

Actual Sample

Five right-handed individuals in the age range 22-27 years participated in the study, four of them female and one male. No exclusion criteria were specified in the analysis plan and no exclusions were necessary.

Differences from pre-data collection methods plan

The verbal instructions to participants were modified after the first participant, to include what they should do if the information tab disappeared; and after the third participant, to include a physical demonstration of a key press. Procedures for dealing with double readings are described above and were used consistently for all 5 participants. Abduction angles for adjacent fingers, measured in the plane containing the metacarpal bones (in the palm) of the same fingers, were added to the list of recorded data before collecting data from the fifth participant; these are used in exploratory analysis only, and not for the confirmatory analysis. Extra rows were created manually in the Excel spreadsheet to accommodate repeated measurements and additional measurements. These measurements were indicated by suffixes to the joint names to allow for filtering during analysis.

Results

The results presented below confirm that the cumulative proportion of variance accounted for by the first two principal components indeed exceeds 80% for each of the five participants. Exploratory plots of the objects in the PC1/PC2 plane are also presented. Participant 1 was able to record 54 of 57 hand postures with 3 lying outside the field of view of the Hololens, while the other 4 participants recorded all 57.

Data preparation

Data preparation following the analysis plan:

###Data Preparation

####Load Relevant Libraries and Functions
library(foreign) # for reading spss formatted data
library(tidyr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stringr) # useful for some string manipulation
library(ggplot2)
library(ggrepel)
library(readxl) # import excel files

####Import data
data_path <- '../participant_data/santello1998_participant_data.xlsx'
d1 = read_excel(data_path, sheet=1)
d1df = data.frame(d1)

#hololens gives data that is easiest to paste in long form

#### Data exclusion / filtering

#### Prepare data for analysis - create columns etc.

handpose_long_unfiltered = select(d1df, -Timestamp_seconds)

#post data collection addendum
#exclude double presses and readings reserved for exploratory analysis
handpose_long = handpose_long_unfiltered %>%
  filter(!str_detect(Joint, "proj")) %>%
  filter(!str_detect(Object, "ignore"))
  #remove repeated observations marked during experiment
  #remove projected abduction angles for confirmatory analysis

handpose_wide_raw = pivot_wider(handpose_long, names_from = Joint, values_from = Angle_degrees)

#average over repeated trials
handpose_means =  handpose_long %>% 
  group_by(ParticipantID, Object, Joint) %>%
  summarise(Mean_angle_degrees = mean(Angle_degrees, na.rm = TRUE)) %>% #this will find mean over trials for given participant and object
  ungroup()
## `summarise()` has grouped output by 'ParticipantID', 'Object'. You can override using the `.groups` argument.
#wide format for all participants
handpose_wide_processed_uncentered = handpose_means %>%
  group_by(ParticipantID)%>%
  pivot_wider(names_from = "Object", values_from = "Mean_angle_degrees")%>%
  ungroup

#wide format for all participants, centered around mean for each joint angle calculated over objects
handpose_wide_processed_all <-
  data.frame(ParticipantID = handpose_wide_processed_uncentered$ParticipantID, Joint = handpose_wide_processed_uncentered$Joint,
  select(handpose_wide_processed_uncentered, -c(ParticipantID, Joint))-rowMeans(select(handpose_wide_processed_uncentered, -c(ParticipantID, Joint)), na.rm = TRUE))

Confirmatory analysis

The single statistical test that justifies the main inference of the paper is whether the amount of variation accounted for by the first two principal components together exceeds 80%.

#split data into participants for individual analysis
datalist = split(handpose_wide_processed_all, f= handpose_wide_processed_all$ParticipantID)

njoints = length(unique(handpose_long$Joint)) #number of joints
nobj = length(unique(handpose_long$Object)) #number of objects
npart = length(unique(handpose_long$ParticipantID)) #number of participants

#outputs to compare
pva_list <- vector(mode = "list", length = npart) #list of vectors containing percent variation explained by each component
coeffs_list <- data.frame(Participant = rep(as.factor(c(1:npart)), each=nobj), Object =  rep(as.character(NA), times=npart*nobj), PC1 = rep(NA, times=npart*nobj), PC2 = rep(NA, times=npart*nobj))  #dataframe containing coefficients of PC1 and PC2 for each object and participant (for plotting)
checkcoeffs = vector(mode = "list", length = npart) #should contain values close to zero if calculations ran correctly 

#do PCA for each participant
for(i in 1:npart){
handpose_wide_processed = as.data.frame(datalist[[i]]) #get data for ith participant in a dataframe

#set row and column names
colnames(handpose_wide_processed) = colnames(handpose_wide_processed_all)
row.names(handpose_wide_processed) <- handpose_wide_processed$Joint

#keep only joint angles in the data, remove columns where data is missing (possible if all trials for an object failed)
handpose_wide_processed <- handpose_wide_processed %>%
  select(-c(ParticipantID, Joint)) %>%
    select_if(~ !any(is.na(.)))

handpose_wide_transposed = t(handpose_wide_processed)
  
#eigenvalues and eigenvectors of covariance matrix
jointscov = cov(handpose_wide_transposed)
eigvals = eigen(jointscov)$values 
eigvecs = eigen(jointscov)$vectors

#create array to store joint vectors for objects expressed in PC basis
df_coeffs <- as.data.frame(handpose_wide_transposed-handpose_wide_transposed)
colnames(df_coeffs) <- c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10","PC11", "PC12", "PC13", "PC14","PC15")

#actual number of objects for which data was captured for this participant  
nobj_i = ncol(handpose_wide_processed)
nobj_i

#fill in coefficients of each PC for each object for this participant
for (objnumber in 1:nobj_i){
  for (pcnumber in 1:njoints) {
  df_coeffs[objnumber,pcnumber] <- sum(eigvecs[,pcnumber]*handpose_wide_processed[,objnumber])
  }
}

start = nobj_i*(i-1)+1
end = nobj_i*i
coeffs_list$Object[start:end] <- rownames(df_coeffs)
coeffs_list$PC1[start:end] <- df_coeffs[,1]
coeffs_list$PC2[start:end] <- df_coeffs[,2]

#check -- this should be close to zero in each entry, it is the object data minus the appropriate linear combination of PCs for that object with the calculated coefficients
df_check <- as.data.frame(handpose_wide_processed)
for (objnumber in 1:nobj_i){
  for (pcnumber in 1:njoints) {
  df_check[,objnumber] <- df_check[,objnumber] - df_coeffs[objnumber, pcnumber]*eigvecs[,pcnumber]
  }
}

#frobenius norm will be close to zero only if all entries are close to zero as they should be
checkcoeffs[[i]] = norm(data.matrix(df_check), "f")

# Proportion of variance explained by each principal component
pva_list[[i]] <- eigvals / sum(eigvals)

}

#plot cumulative proportion of variance explained by PCs for all participants
xValue <- 1:njoints
dfplot <- data.frame(Participant = rep(as.factor(c(1:npart)), each=njoints), xValue, c(cumsum(pva_list[[1]]), cumsum(pva_list[[2]]), cumsum(pva_list[[3]]), cumsum(pva_list[[4]]), cumsum(pva_list[[5]])))
colnames(dfplot)<-c("Participant", "PC", "Percent_variance_explained")

vp <- ggplot(dfplot, aes(x = PC, y= Percent_variance_explained, group = Participant, color=Participant)) +
  geom_line() +ggtitle("Cumulative Proportion of Variance Explained")+
  geom_point() + ylab("Cumulative Proportion of Variance Explained") + scale_x_discrete(name="Principal Component", limits = colnames(df_coeffs)) + theme_minimal()+theme(legend.position="top") # + ylim(0, 1)

vp

From the above, the key test is successful, as the cumulative proportion of variance at the second principal component is above 0.8 for every one of the participants. For specific comparison with the original paper, the percent variances accounted for by the first two principal components are presented below (data from original study are from Table 2 of the paper). Of course, the participant IDs are not relevant when comparing across experiments. These data are presented here only to show the details of the percent variation explained.

pva_table <- data.frame(Participant = c(1:npart), OriginalPC1 = c(1:5), OriginalPC2 =  c(1:5), ReplicationPC1 = c(1:5), ReplicationPC2 =  c(1:5), Totalorig =c(1:5),  Totalrep=  c(1:5)) 
pva_table$OriginalPC1 = c(52.9, 49.5, 74.8, 79.3, 62.9)
pva_table$OriginalPC2 = c(24.7, 37.6, 13.0, 10.0, 17.2)
pva_table$Totalorig =pva_table$OriginalPC1+pva_table$OriginalPC2
pva_table$ReplicationPC1 = 100*c(pva_list[[1]][1], pva_list[[2]][1], pva_list[[3]][1], pva_list[[4]][1], pva_list[[5]][1])
pva_table$ReplicationPC2 = 100*c(pva_list[[1]][2], pva_list[[2]][2], pva_list[[3]][2], pva_list[[4]][2], pva_list[[5]][2])
pva_table$Totalrep = pva_table$ReplicationPC1 + pva_table$ReplicationPC2
colnames(pva_table) = c("Participant", "PC1 (original)", "PC2 (original)", "PC1 (replication)", "PC2 (replication)", "PC1+PC2 (original)", "PC1+PC2 (replication)")
pva_table
##   Participant PC1 (original) PC2 (original) PC1 (replication) PC2 (replication)
## 1           1           52.9           24.7          66.12685         14.013787
## 2           2           49.5           37.6          70.61454         13.990449
## 3           3           74.8           13.0          75.33583          8.596148
## 4           4           79.3           10.0          81.35588          8.729056
## 5           5           62.9           17.2          80.98138          6.798513
##   PC1+PC2 (original) PC1+PC2 (replication)
## 1               77.6              80.14064
## 2               87.1              84.60499
## 3               87.8              83.93198
## 4               89.3              90.08494
## 5               80.1              87.77989
#knitr::opts_chunk$set(fig.width = 20) 

#plot objects in the PC1/PC2 space for all participants
theme<-theme(panel.background = element_blank(),panel.border=element_rect(fill=NA),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),strip.background=element_blank(),axis.text.x=element_text(colour="black"),axis.text.y=element_text(colour="black"),axis.ticks=element_line(colour="black"),plot.margin=unit(c(1,1,1,1),"line"))

coeffs_list<- na.omit(coeffs_list)
                                                                                                               
p<-ggplot(coeffs_list,aes(x=PC1, y=PC2, group = Participant))
p<-p+geom_point()+ geom_text_repel(size=5, label=coeffs_list$Object, max.overlaps = Inf)+theme+facet_wrap(Participant~., ncol=1, labeller = label_both) + scale_x_reverse() + scale_y_reverse() +coord_fixed(ratio = 1)
p

#edited output figure size of the code chunk (thanks Xi Jia!) and font size of ggrepel post data collection

#AFTER DATA COLLECTION, some of the axes of PC plots may need to be flipped just to make visual comparison easier (sign of eigenvector doesn't matter to the analysis)

The plot of objects in the PC1/PC2 plane for a single participant is given in the original paper as Figure 7 of that paper. It is reproduced below for comparison. While the scale of coefficients on the x-axis is very different in the replication, some patterns seem to be common to the two studies. For example, at least the plots for participants 3, 4, and 5 seem to have a bilinear distribution. Interestingly, the authors of the original study also found such a pattern in exactly 3 of their 5 subjects, including the one whose data is the source of Figure 7.

Comparing the locations of specific objects in the plot from the original paper with Participant 5 from the replication, the relative positions of the CD, screwdriver, coffee mug, baseball, egg, computer mouse, wrench, milk carton, notebook, tennis racket, dictionary, doorknob, frying pan, knob of a lid, cherry, dinner plate and several others seem to match quite well between the two plots, despite the mismatch in numerical values (which is likely due to some scaling of unit eigenvectors by the original authors which was not done in the replication). At least, objects that are together on the right or left side of one plot seem to be on the same side in the other plot as well. Examples of objects whose locations do not match are the Chinese teacup and turtle, both less universally-encountered objects to which participants in the replication expressed bewilderment. Some common objects such as the zipper and chalk are also out of alignment. Still, the existence of any correspondence at all in this context is surprising, particularly since there were differences across participants in how the prompts were interpreted, and at least two participants mentioned informally that they may have imagined different grasps in different trials of the same object. If indeed such a grouping (separation by the sign of the first principal component) is significant, it would be very interesting to examine what it represents. As pointed out by the authors of Santello et al., it does not appear to be related to grasp taxonomies.

#![Plot of objects in the PC1/PC2 plane from the original paper (Fig. 7)](fig7_santello1998.jpg)
knitr::include_graphics("fig7_santello1998.jpg")

Exploratory analyses

One of the exploratory analyses is shown above as it was planned prior to the data collection. This section includes two other exploratory analyses: the first, replacement of the angles between the proximal phalanges with the angle between their projections on the palm (potentially a closer measure of abduction angle to that used in the original study), and the second, a repetition of the confirmatory test after scaling of joint angles by their standard deviations. This scaling step is common in PCA, however, it is not mentioned in the original paper despite centering about the mean being explicitly mentioned. For that reason, it was omitted in the confirmatory analysis.

#Using projections of fingers to find abduction angles in the plane of the palm -- for Participant 5 alone

handpose_long_proj = handpose_long_unfiltered %>%
    filter(!str_detect(Object, "ignore"))%>%
  #remove repeated observations marked during experiment
  filter(ParticipantID == "5" & (!str_detect(Joint, "ABD") | (str_detect(Joint, "proj") | str_detect(Joint, "Thumb"))))
  #select data belonging to participant 5, which is either not an abduction angle, or if it is, then is either of the thumb or a projected angle

handpose_wide_raw_proj = pivot_wider(handpose_long_proj, names_from = Joint, values_from = Angle_degrees)

#average over repeated trials
handpose_means_proj =  handpose_long_proj %>% 
  select(-c(ParticipantID))%>%
  group_by(Object, Joint) %>%
  summarise(Mean_angle_degrees = mean(Angle_degrees, na.rm = TRUE)) %>% #this will find mean over trials for given participant and object
  ungroup()
## `summarise()` has grouped output by 'Object'. You can override using the `.groups` argument.
#wide format 
handpose_wide_processed_uncentered_proj = handpose_means_proj %>%
  pivot_wider(names_from = "Object", values_from = "Mean_angle_degrees")%>%
  ungroup

#wide format for all participants, centered around mean for each joint angle calculated over objects
handpose_wide_processed_proj <-
  data.frame(Joint = handpose_wide_processed_uncentered_proj$Joint,
  select(handpose_wide_processed_uncentered_proj, -c(Joint))-rowMeans(select(handpose_wide_processed_uncentered_proj, -c(Joint)), na.rm = TRUE))

njoints5 = length(unique(handpose_long_proj$Joint)) #number of joints
nobj5 = length(unique(handpose_long_proj$Object)) #number of objects

#set row and column names
row.names(handpose_wide_processed_proj) <- handpose_wide_processed_proj$Joint

#keep only joint angles in the data, remove columns where data is missing (possible if all trials for an object failed)
handpose_wide_processed_proj <- handpose_wide_processed_proj %>%
  select(-c(Joint)) %>%
    select_if(~ !any(is.na(.)))
handpose_wide_transposed_proj = t(handpose_wide_processed_proj)
  
#eigenvalues and eigenvectors of covariance matrix
jointscov = cov(handpose_wide_transposed_proj)
eigvals = eigen(jointscov)$values 
eigvecs = eigen(jointscov)$vectors

#create array to store joint vectors for objects expressed in PC basis
df_coeffs5 <- as.data.frame(handpose_wide_transposed_proj-handpose_wide_transposed_proj)
colnames(df_coeffs5) <- c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10","PC11", "PC12", "PC13", "PC14","PC15")

#actual number of objects for which data was captured for this participant  
nobj5 = ncol(handpose_wide_processed_proj)

#fill in coefficients of each PC for each object for this participant
for (objnumber in 1:nobj5){
  for (pcnumber in 1:njoints) {
  df_coeffs5[objnumber,pcnumber] <- sum(eigvecs[,pcnumber]*handpose_wide_processed_proj[,objnumber])
  }
}

coeffs_list5 <- data.frame(Object =  rep(as.character(NA), times=nobj5), PC1 = rep(NA, times=nobj5), PC2 = rep(NA, times=nobj5))  #dataframe containing coefficients of PC1 and PC2 for each object (for plotting)
start = 1
end = nobj5
coeffs_list5$Object[start:end] <- rownames(df_coeffs5)
coeffs_list5$PC1[start:end] <- df_coeffs5[,1]
coeffs_list5$PC2[start:end] <- df_coeffs5[,2]

#check -- this should be close to zero in each entry, it is the object data minus the appropriate linear combination of PCs for that object with the calculated coefficients
df_check5 <- as.data.frame(handpose_wide_processed_proj)
for (objnumber in 1:nobj5){
  for (pcnumber in 1:njoints) {
  df_check5[,objnumber] <- df_check5[,objnumber] - df_coeffs5[objnumber, pcnumber]*eigvecs[,pcnumber]
  }
}


#frobenius norm will be close to zero only if all entries are close to zero as they should be
norm(data.matrix(df_check5), "f")
## [1] 1.203422e-12
# Proportion of variance explained by each principal component
pva_list5 <- eigvals / sum(eigvals)

#plot cumulative proportion of variance explained by PCs
dfplot5 <- data.frame(Analysis = rep(as.factor(c(1:2)), each = njoints5), xval = rep(xValue, times=2), pva = c(cumsum(pva_list[[5]]), cumsum(pva_list5)))
colnames(dfplot5)<-c("Analysis", "PC", "PVE")

vp5 <- ggplot(dfplot5, aes(x = PC, y= PVE, group = Analysis, color = Analysis)) +
  geom_line() +ggtitle("Cumulative Proportion of Variance Explained for Participant 5")+
  geom_point() + ylab("Cumulative Proportion of Variance Explained") + scale_x_discrete(name="Principal Component", limits = colnames(df_coeffs5)) + theme_minimal()

vp5

#plot objects in the PC1/PC2 space for participant 5
coeffs_list5<- na.omit(coeffs_list5)

expplot1 <- data.frame(Analysis = rep(as.factor(c(1:2)), each=nobj5), Object = rep(coeffs_list5$Object, times=2), PC1 = c(filter(coeffs_list, coeffs_list$Participant=="5")$PC1, (-1)*coeffs_list5$PC1), PC2 = c(filter(coeffs_list, coeffs_list$Participant=="5")$PC2, coeffs_list5$PC2))

p512<-ggplot(expplot1, aes(x=PC1, y=PC2, group = Analysis, color = Analysis))
p512<-p512+geom_point()+ geom_text_repel(size=2, label=expplot1$Object, max.overlaps = Inf)+theme + scale_x_reverse() + scale_y_reverse() +coord_fixed(ratio = 1) # + facet_wrap(Analysis~., ncol=2, labeller = label_both) + 
p512

With this method of measurement, the plots of cumulative proportion of variance and objects in the PC1/PC2 plane both shift slightly, however, the first two principal components continue to account for more than 0.8 of the variance.

#with scaling by standard deviation

handpose_wide_transposed = handpose_wide_transposed %>%
  apply(2,function(x){x/sd(x)}) #scaling step

  
#eigenvalues and eigenvectors of covariance matrix
jointscov = cov(handpose_wide_transposed)
eigvals = eigen(jointscov)$values 
eigvecs = eigen(jointscov)$vectors

#create array to store joint vectors for objects expressed in PC basis
df_coeffs5 <- as.data.frame(handpose_wide_transposed-handpose_wide_transposed)
colnames(df_coeffs5) <- c("PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10","PC11", "PC12", "PC13", "PC14","PC15")

#actual number of objects for which data was captured for this participant  
nobj5 = ncol(handpose_wide_processed)

#fill in coefficients of each PC for each object for this participant
for (objnumber in 1:nobj5){
  for (pcnumber in 1:njoints) {
  df_coeffs5[objnumber,pcnumber] <- sum(eigvecs[,pcnumber]*t(handpose_wide_transposed)[,objnumber])
  }
}

coeffs_list5 <- data.frame(Object =  rep(as.character(NA), times=nobj), PC1 = rep(NA, times=nobj), PC2 = rep(NA, times=nobj))  #dataframe containing coefficients of PC1 and PC2 for each object (for plotting)
start = 1
end = nobj5
coeffs_list5$Object[start:end] <- rownames(df_coeffs5)
coeffs_list5$PC1[start:end] <- df_coeffs5[,1]
coeffs_list5$PC2[start:end] <- df_coeffs5[,2]

#check -- this should be close to zero in each entry, it is the object data minus the appropriate linear combination of PCs for that object with the calculated coefficients
df_check5 <- as.data.frame(t(handpose_wide_transposed))
for (objnumber in 1:nobj5){
  for (pcnumber in 1:njoints) {
  df_check5[,objnumber] <- df_check5[,objnumber] - df_coeffs5[objnumber, pcnumber]*eigvecs[,pcnumber]
  }
}


#frobenius norm will be close to zero only if all entries are close to zero as they should be
norm(data.matrix(df_check5), "f")
## [1] 2.292521e-14
# Proportion of variance explained by each principal component
pva_list5 <- eigvals / sum(eigvals)

#plot cumulative proportion of variance explained by PCs 
dfplot5 <- data.frame(Analysis = rep(as.factor(c(1:2)), each = njoints), xval = rep(xValue, times=2), pva = c(cumsum(pva_list[[5]]), cumsum(pva_list5)))
colnames(dfplot5)<-c("Analysis", "PC", "PVE")

vp5 <- ggplot(dfplot5, aes(x = PC, y= PVE, group = Analysis, color = Analysis)) +
  geom_line() +ggtitle("Cumulative Proportion of Variance Explained for Participant 5")+
  geom_point() + ylab("Cumulative Proportion of Variance Explained") + scale_x_discrete(name="Principal Component", limits = colnames(df_coeffs5)) + theme_minimal()

vp5

#plot objects in the PC1/PC2 space for participant 5
coeffs_list5<- na.omit(coeffs_list5)

expplot1 <- data.frame(Analysis = rep(as.factor(c(1:2)), each=nobj5), Object = rep(coeffs_list5$Object, times=2), PC1 = c(filter(coeffs_list, coeffs_list$Participant=="5")$PC1, coeffs_list5$PC1), PC2 = c(filter(coeffs_list, coeffs_list$Participant=="5")$PC2, coeffs_list5$PC2))

p512<-ggplot(expplot1, aes(x=PC1, y=PC2, group = Analysis, color = Analysis))
p512<-p512+geom_point()+ geom_text_repel(size=3, label=expplot1$Object, max.overlaps = Inf)+ scale_x_reverse() + scale_y_reverse() + facet_wrap(Analysis~., ncol=2, labeller = label_both, scales="free") 
p512

Scaling the readings for each joint by their standard deviation reduces the amount of variance explained by the first two principal components, and the cumulative proportion of variance does not exceed 0.8 until as many as four principal components are included. The plot of coefficients shrinks by almost an order of magnitude along both PC1 and PC2 axes.

It is worth mentioning here that PC1 being a principal component implies that -PC1 is also the same principal component, therefore, axes of the PC1/PC2 plots have been inverted on several occasions within this document, and at times, the values been plotted have been all negated at once to make comparison with another plot easier.

Discussion

The primary result, that the first two principal components account for more than 80% of the variance in joint angles in all participants, was successfully replicated according to the measurement method defined at the outset. The exploratory analysis suggests that changing the way in which abduction angles are estimated will not affect this result.

It stands to reason that an inference about the effective degrees-of-freedom of the hand should not be affected by the specific parameters used to represent hand pose. It is encouraging to note that, as expected, the result is robust to a change of parameters, a change of measurement device, and the passage of 23 years. Certain objects in the list have lost their relevance since the original paper was published (e.g, compact disc); some were functionally unfamiliar to some replication participants (ashtray, Chinese tea cup); some were imagined in various shapes by various participants (e.g. drawer handle) and even across trials by the same participant, some were ``used" in different ways (e.g., lifting or patting a turtle), some typically require bimanual operation and were consequently treated differently by different participants (e.g., fishing rod). Yet the diversity of the list remains adequate to test the hypothesis and arrive at the same conclusion. Moreover, for many objects, the arrangements of objects along the PC1/PC2 plane that are reported in the original paper reappear in the replication.

It is exciting to note this last fact, that objects seem to group together on either side of a vertical axis (or, along two oblique lines) in a very similar manner across participants, including when compared to the participant whose data is plotted in the original paper. If these groupings are unrelated to grasp taxonomy, it must at least have meaning with respect to the handshapes represented by the two principal components. The original paper includes a discussion of the continuum between poses represented by moving on the PC1/PC2 plane; a good extension of this replication project would be to compare the principal components derived to those in the original paper, in terms of the handshape that they represent, and to examine – by visualizing hand postures – why the remaining two subjects’ postures do not seem to be arranged in the same way. Another extension would be to complete the reproduction of other results provided in the paper, such as the information theoretic analysis that led the authors to assess 5 degrees of freedom as being sufficient to capture fine adjustments in hand posture.