As mentioned in the data exploration document, eccentricity did not seem normally distributed, and could not be effectively transformed to resemble a more Gaussian distribution. In this document, I attempted to model eccentricity using different non-linear functions such that I could extract the relevant coefficients/parameters for those function for model-based clustering.
There are different non-linear functions we can use to fit the data. I fit the data using 4 different types of non-linear functions to see which one might be the most appropriate to use
# ----------- Load libraries ------------
library(plyr)
library(openxlsx)
# ------------ Load data -------------
df<-read.csv("VR_Basic Scene Viewing_Summary Per Trial_Updated 11Nov2020.csv")
suj<-unique(df$subject)
ecc<-NULL
head.coef<-NULL
sacc.coef<-NULL
# definining non-linear functions
f.exp <- function(X,a,b) {a * (X^(b - 1) * b)} # exponential function
f.pwr <- function(X,a,k) {a * (exp(k * X) * k)} # power function
for (s in suj){
# subset data frame for each participant
sub<-subset(df, subject==s)
# sort eccentricities from smallest to largest and create data frame
frame<-data.frame(subject = s,
len = 1:nrow(sub),
head.ecc = sort(sub$head.ecc),
sacc.ecc = sort(sub$sacc.ecc)
)
# bind data frame
ecc<-rbind(ecc, frame)
# ---------- HEAD ------------
# EXPONENTIAL FUNCTION
# get the starting position of a and b using a linear model
fexp0 <- nls(log(head.ecc) ~ log(f.exp(len, a, b)), frame, start = c(a = 1, b = 1))
# construct exponential model and extracting coefficients
mod.exp<-nls(head.ecc ~ f.exp(len, a, b), frame, start = coef(fexp0))
mod.exp.sum<-summary(mod.exp)
# POWER FUNCTION
# get the starting position
fpwr0 <- nls(log(head.ecc) ~ log(f.pwr(len, a, k)), frame, start = c(a = 1, k = 1))
# construct model and extracting coefficients
mod.pwr<-nls(head.ecc ~ f.pwr(len, a, k), data=frame, start = coef(fpwr0))
mod.pwr.sum<-summary(mod.pwr)
# ASYMPTOTIC REGRESSION
mod.asym<-nls(head.ecc ~ SSasymp(len, Asym, R0, lrc), data=frame)
mod.asym.sum<-summary(mod.asym)
# RECTANGULAR HYPERBOLA
mod.mm<-nls(head.ecc ~ SSmicmen(len, Vm, K), data=frame)
mod.mm.sum<-summary(mod.mm)
h.coef<-data.frame(subject = s,
# exponential function coef
exp.int = mod.exp.sum$coefficients[1],
exp.rate = mod.exp.sum$coefficients[2],
# power function coef
pwr.int = mod.pwr.sum$coefficients[1],
pwr.rate = mod.pwr.sum$coefficients[2],
# asymptotic regression
asym.max = mod.asym.sum$coefficients[1],
asym.int = mod.asym.sum$coefficients[2],
asym.lograte = mod.asym.sum$coefficients[3],
# Rectangular hyperbola
mm.max = mod.mm.sum$coefficients[1],
mm.maxhalf = mod.mm.sum$coefficients[2],
# fit between model and data
exp.cor = cor(frame$head.ecc, predict(mod.exp)),
pwr.cor = cor(frame$head.ecc, predict(mod.pwr)),
asym.cor = cor(frame$head.ecc, predict(mod.asym)),
mm.cor = cor(frame$head.ecc, predict(mod.mm))
)
head.coef<-rbind(head.coef, h.coef)
# ---------- SACCADE ------------
# EXPONENTIAL FUNCTION
# get the starting position of a and b using a linear model
fexp0 <- nls(log(sacc.ecc) ~ log(f.exp(len, a, b)), data=frame, start = c(a = 1, b = 1))
# construct exponential model and extracting coefficients
mod.exp<-nls(sacc.ecc ~ f.exp(len, a, b), data=frame, start = coef(fexp0))
mod.exp.sum<-summary(mod.exp)
# POWER FUNCTION
# get the starting position
fpwr0 <- nls(log(sacc.ecc) ~ log(f.pwr(len, a, k)), data=frame, start = c(a = 1, k = 1))
# construct model and extracting coefficients
mod.pwr<-nls(sacc.ecc ~ f.pwr(len, a, k), data=frame, start = coef(fpwr0))
mod.pwr.sum<-summary(mod.pwr)
# ASYMPTOTIC REGRESSION (using self-starting function)
mod.asym<-nls(sacc.ecc ~ SSasymp(len, Asym, R0, lrc), data=frame)
mod.asym.sum<-summary(mod.asym)
# RECTANGULAR HYPERBOLA (using self-starting function)
mod.mm<-nls(sacc.ecc ~ SSmicmen(len, Vm, K), data=frame)
mod.mm.sum<-summary(mod.mm)
s.coef<-data.frame(subject = s,
# exponential function coef
exp.int = mod.exp.sum$coefficients[1],
exp.rate = mod.exp.sum$coefficients[2],
# power function coef
pwr.int = mod.pwr.sum$coefficients[1],
pwr.rate = mod.pwr.sum$coefficients[2],
# asymptotic regression
asym.max = mod.asym.sum$coefficients[1],
asym.int = mod.asym.sum$coefficients[2],
asym.lograte = mod.asym.sum$coefficients[3],
# Rectangular hyperbola
mm.max = mod.mm.sum$coefficients[1],
mm.maxhalf = mod.mm.sum$coefficients[2],
# fit between model and data
exp.cor = cor(frame$sacc.ecc, predict(mod.exp)),
pwr.cor = cor(frame$sacc.ecc, predict(mod.pwr)),
asym.cor = cor(frame$sacc.ecc, predict(mod.asym)),
mm.cor = cor(frame$sacc.ecc, predict(mod.mm))
)
sacc.coef<-rbind(sacc.coef, s.coef)
}
# datasets <- list("Saccade Paramters" = sacc.coef, "Head Movement Parameters" = head.coef)
# write.xlsx(datasets, file = "Eccentricity Parameters from Non-Linear Models.xlsx", row.names=F)
# ----------- creating summary tables ----------------
# Mean correlation between actual and model value
head.coef.sum<-ddply(head.coef, .(), summarise,
exp = mean(exp.cor),
pwr = mean(pwr.cor),
asym = mean(asym.cor),
mm = mean(mm.cor))
# Mean correlation between actual and model value
sacc.coef.sum<-ddply(sacc.coef, .(), summarise,
exp = mean(exp.cor),
pwr = mean(pwr.cor),
asym = mean(asym.cor),
mm = mean(mm.cor))Below are tables examining how well the predicted values fit with actual values using the correlation function. As shown, the asymptotic regression model might fit the data best for both eye and head movement eccentricity.
## exp pwr asym mm
## 1 0.9511513 0.8117374 0.9754485 0.9218623## exp pwr asym mm
## 1 0.939488 0.7794861 0.9746239 0.9314913