Last active
March 25, 2018 08:36
-
-
Save kagaya/cbafe0a332e6d6766168 to your computer and use it in GitHub Desktop.
functions for statistical analysis and data visualizations in the paper — Motor control of ultrafast, ballistic movements (by K. Kagaya and S. N. Patek)
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| ## functions for statistical analysis and data visualizations in the paper, | |
| ## The paper's title: Motor control of ultrafast, ballistic movements (by K. Kagaya and S. N. Patek) | |
| ## This script is written by Katsushi Kagaya ([email protected]) | |
| ## Gist URL: https://gist.github.com/kagaya/cbafe0a332e6d6766168 | |
| ## to load libraries | |
| library(ggplot2) # for plotting data and model prediction | |
| library(mgcv) # for generalized addive modeling (GAM fitting) | |
| library(plyr) # for data manipulation | |
| library(nlme) # for linear mixed effects modeling | |
| generate.data.for.model.selection <- | |
| function(ang = 20, time.lim = c(-10, 0), ...) { | |
| ## to generate data set (stat.df) for statistical modeling | |
| ## angular positon for measuring instantaneous velocity | |
| ## limits for the time window selecting spikes, the unit is second | |
| ## give data sets, like, "nb18, nb30, nb55, nb500, nb504, nb520" | |
| arg.len <- length(list(...)) | |
| df <- data.frame() | |
| for (i in 1:arg.len) { | |
| if (sum(is.element(names(list(...)[[i]]), "strike")) == 1) { | |
| raw <- list(...)[[i]]$strike$raw | |
| r <- list(...)[[i]]$strike$output.propodus.rotation | |
| r.meralV <- list(...)[[i]]$strike$output.meralV.rotation | |
| r.in <- list(...)[[i]]$contraction$input.rotation | |
| spike.df <- list(...)[[i]]$spike.df | |
| fps <- list(...)[[i]]$fps.strike | |
| fps.c <- list(...)[[i]]$fps.contraction | |
| mmperpix <- list(...)[[i]]$mmperpix | |
| ## extract meral-V kinematics and store as data frame | |
| r.in2 <- data.frame() | |
| for (ii in 1:length(r.in)) { | |
| x <- -seq(0, 1 / fps.c * (length(r.in[[ii]]) - 1), 1 / fps.c) | |
| x <- x[length(x):1] | |
| r.in2 <- rbind(r.in2, | |
| data.frame( | |
| strike.num = ii, time = x, rotation = r.in[[ii]] )) | |
| } | |
| on <- ddply(spike.df, .(strike.num), subset, time == min(time)) | |
| rot.sub <- data.frame() | |
| init.rot <- c() | |
| for (ii in on$strike.num) { | |
| sub.df <- subset(r.in2, | |
| strike.num == ii & | |
| time >= on$time[ii]) # rotation after extensor onset | |
| init.rot <- c(init.rot, sub.df$rotation[1]) | |
| sub.df$rotation <- sub.df$rotation - sub.df$rotation[1] | |
| rot.sub <- rbind(rot.sub, sub.df) | |
| } | |
| r.in2.df <- data.frame() | |
| for (ii in 1:length(init.rot)) { | |
| sub.rot.df <- subset(r.in2, strike.num == ii) | |
| sub.rot.df$rotation <- | |
| sub.rot.df$rotation - init.rot[ii] | |
| r.in2.df <- rbind(r.in2.df, sub.rot.df) | |
| } | |
| ## make a data frame for the storage of kinematic data | |
| kinem.df <- data.frame() | |
| for (j in 1:length(r)) { | |
| if (length(r[[j]] > 0)) { | |
| # because there is missing data in propodus rotation | |
| ## output propodus rotation | |
| time.temp <- seq(0, 1 / fps * (length(r[[j]]) - 1), 1 / fps) | |
| rot <- data.frame(strike.num = j, time = time.temp, rotation = r[[j]]) | |
| fit <- gam(rotation ~ s(time, k = 5), data = rot) | |
| time.grid <- seq(0, max(rot$time), 1/fps/10) | |
| pred <- predict(fit, newdata = data.frame(time = time.grid), type = "response") | |
| rot2 <- subset(data.frame(time = time.grid, rotation = pred), rotation <= ang) | |
| ## time at certain rotation degrees set by ang | |
| time.at.limit <- rot2$time[length(rot2$time)] | |
| vel.temp <- diff(rot2$rotation) / (1/fps/10) | |
| vel2 <- data.frame(time = rot2$time, velocity = c(vel.temp[1], vel.temp)) | |
| ## maximum angular velocity | |
| max.vel <- max(vel.temp) | |
| acc.temp <- diff(vel2$velocity) / (1/fps/10) | |
| acc <- data.frame(time = rot2$time, acceleration = c(acc.temp[1], acc.temp)) | |
| ## maximum acceleration | |
| max.acc <- max(acc$acceleration) | |
| ## acceleration at limit | |
| acc.at.limit <- acc.temp[length(acc.temp)] | |
| ## velocity at limit ( 1/fps/10) | |
| vel.at.limit <- vel.temp[length(vel.temp)] | |
| ## time at maximum acceleration | |
| time.at.max.acc <- subset(acc, acceleration == max.acc)$time[1] | |
| ## linear velocity | |
| point.x <- subset(raw[[j]], track.id == 1)$point.x | |
| point.y <- subset(raw[[j]], track.id == 1)$point.y | |
| linear.distance <- cumsum(sqrt(diff(point.x) ^ 2 + diff(point.y) ^ 2) * mmperpix) # mm | |
| time.temp <- seq(0, 1 / fps * (length(point.x) - 2), 1 / fps) | |
| dist <- data.frame(time = time.temp, linear.distance = linear.distance) | |
| fit.dist <- gam(linear.distance ~ s(time, k = 5), data = dist) | |
| time.grid <- seq(0, max(dist$time), 1/fps/10) # for visualization, 1/10 time step | |
| pred.dist <- predict(fit.dist, newdata = data.frame(time = time.grid), type = "response") | |
| dist2 <- | |
| subset( | |
| data.frame(time = time.grid, linear.distance = pred.dist), | |
| time <= time.at.limit | |
| ) | |
| linear.vel.temp <- diff(dist2$linear.distance) / (1/fps/10) / 1000 # m/s | |
| ## linear velocity at limit (1/fps/10) | |
| linear.vel.at.limit <- linear.vel.temp[length(linear.vel.temp)] | |
| linear.vel.all <- | |
| data.frame( | |
| time = dist2$time, linear.velocity = c(linear.vel.temp[1], linear.vel.temp) | |
| ) | |
| ## linear acceleration | |
| linear.acc.temp <- diff(linear.vel.all$linear.velocity) / (1/fps/10) # m/s^2 | |
| ## ## maximum linear acc | |
| max.linear.acc <- max(linear.acc.temp) | |
| ## ## linear acceleration at limit | |
| linear.acc.at.limit <- linear.acc.temp[length(linear.acc.temp)] # m/s^2 | |
| ## maximum linear velocity | |
| max.linear.vel <- max(linear.vel.all$linear.velocity) # m/s | |
| ## duration to maximum linear velocity | |
| time.at.max.linear.vel <- | |
| linear.vel.all$time[linear.vel.all$linear.velocity == max.linear.vel] * 10 ^ 3 # ms | |
| ## no missing data in input rotation, but reduce the data to match | |
| ## the number of data to propodus rotation | |
| input.rotation <- abs(r.in[[j]]) | |
| input.rotation <- input.rotation[length(input.rotation)] | |
| ## input rotation offset by the value when the extensor spike on | |
| input.rotation2 <- abs(min(subset(r.in2.df, strike.num == j)$rotation)) | |
| ## input meral-V rotation, velocity, acceleration | |
| mr.in <- -r.in[[j]] | |
| time.temp <- | |
| seq(0, 1 / fps.c * (length(mr.in) - 1), 1/fps.c) | |
| mr.in.d <- data.frame(time = time.temp, rotation = mr.in) | |
| fit <- gam(rotation ~ s(time, k = 5), data = mr.in.d) | |
| time.grid <- seq(0, max(mr.in.d$time), 1/fps.c/10) | |
| pred <- predict(fit, newdata = data.frame(time = time.grid), type = "response") | |
| mr.in.d2 <- data.frame(time = time.grid, rotation = pred) | |
| vel.temp <- diff(mr.in.d2$rotation) / (1/fps.c/10) | |
| rotational.vel.m.in <- max(vel.temp) | |
| vel.in.d <- data.frame(time = time.grid, | |
| velocity = c(vel.temp[1], vel.temp)) | |
| acc.temp <- diff(vel.in.d$velocity) / (1/fps.c/10) | |
| max.acc.m.in <- max(acc.temp) | |
| ## output meral-V rotation | |
| xx <- seq(0, 1 / fps * (length(r.meralV[[j]]) - 1), 1/fps) | |
| mrd <- data.frame(strike.num = j, time = xx, rotation = r.meralV[[j]]) | |
| fit2 <- gam(rotation ~ s(time, k = 5), data = mrd) | |
| time.grid2 <- seq(0, max(mrd$time), 1/fps/10) | |
| pred2 <- predict(fit2, newdata = data.frame(time = time.grid2), type = "response") | |
| mrd2 <- | |
| subset(data.frame(time = time.grid2, rotation = pred2), | |
| time <= time.at.limit) # time at a certain degree set by ang | |
| vel.meralV.temp <- diff(mrd2$rotation) / (1/fps/10) | |
| # maximum velocity of meralV of unloading | |
| vel.meralV <- max(vel.meralV.temp) | |
| ## velocity of meral-V can not be negative value, | |
| ## but there is in nb500 a negative value (strike 12). | |
| ## It is likely to be generated by an artifuactual appendage movement. | |
| if (vel.meralV < 0) { | |
| vel.meralV <- NA | |
| } | |
| vel2.meralV <- data.frame( | |
| time = mrd2$time, | |
| velocity = c(vel.meralV.temp[1], vel.meralV.temp) | |
| ) | |
| acc.meralV.temp <- | |
| diff(vel2.meralV$velocity) / (1/fps/10) | |
| acc.meralV <- data.frame( | |
| time = mrd2$time, | |
| acceleration = c(acc.meralV.temp[1], acc.meralV.temp) | |
| ) | |
| max.acc.meralV <- max(acc.meralV$acceleration) | |
| time.at.max.acc.meralV <- | |
| subset(acc.meralV, acceleration == max.acc.meralV)$time[1] | |
| sub.kinem.df <- data.frame( | |
| strike.num = j, | |
| ## meralV kinematics | |
| input.rotation = input.rotation, | |
| input.rotation2 = input.rotation2, | |
| rotational.vel.m.in = rotational.vel.m.in, | |
| max.acc.m.in = max.acc.m.in, | |
| velocity.meralV = as.numeric(vel.meralV), | |
| max.acc.meralV = max.acc.meralV, | |
| time.at.max.acc.meralV = time.at.max.acc.meralV, | |
| ## strike duration | |
| time.at.limit = time.at.limit, | |
| time.at.max.linear.vel = time.at.max.linear.vel, # m/s | |
| time.at.max.acc = time.at.max.acc, | |
| ## strike linear kinematics | |
| max.linear.vel = max.linear.vel, | |
| linear.vel.at.limit = linear.vel.at.limit, | |
| max.linear.acc = max.linear.acc, | |
| linear.acc.at.limit = linear.acc.at.limit, | |
| ## strike angular kinematics | |
| max.vel = max.vel, | |
| vel.at.limit = vel.at.limit, | |
| max.acc = max.acc, | |
| acc.at.limit = acc.at.limit | |
| ) | |
| kinem.df <- rbind(kinem.df, sub.kinem.df) | |
| } | |
| } | |
| } | |
| ## generating dataset for EMG spikes | |
| spike.df <- list(...)[[i]]$spike.df | |
| spike.df <- | |
| spike.df[spike.df$time >= time.lim[1] & | |
| spike.df$time <= time.lim[2],] | |
| spike.num <- as.numeric(table(spike.df$strike.num)) | |
| ## no kinematic data! ## | |
| if (sum(is.element(names(list(...)[[i]]), "strike")) == 0) { | |
| colnames.kinem.df <- colnames(kinem.df) | |
| kinem.arr <- array(NA, c(1, length(colnames.kinem.df))) | |
| kinem.df <- as.data.frame(kinem.arr) | |
| colnames(kinem.df) <- colnames.kinem.df | |
| } | |
| ## apply the strike IDs whose number is smaller than the other | |
| ## and save the IDs into strike.id | |
| if (dim(kinem.df)[1] == 1) { | |
| ## when no kinematic data | |
| spike.num <- spike.num[unique(spike.df$strike.num)] | |
| strike.id <- unique(spike.df$strike.num) | |
| } else if (length(kinem.df$strike.num) < length(spike.df$strike.num)) { | |
| ## # strikes is smaller than # EMG | |
| spike.num <- spike.num[kinem.df$strike.num] | |
| strike.id <- kinem.df$strike.num | |
| } else { | |
| ## other case | |
| spike.num <- spike.num[spike.df$strike.num] | |
| strike.id <- spike.df$strike.num | |
| cat("## other case ##\n") | |
| } | |
| dmax <- | |
| ddply(spike.df, .(strike.num), subset, time == max(time)) | |
| dmin <- | |
| ddply(spike.df, .(strike.num), subset, time == min(time)) | |
| ddif <- dmax - dmin | |
| ddif$strike.num <- dmax$strike.num | |
| duration <- | |
| ddif$time # duration does not include silent phase <=> co-activation | |
| duration <- duration[strike.id] | |
| ddif2 <- dmin | |
| ddif2$time <- abs(ddif2$time) | |
| duration2 <- ddif2$time # duration2 includes silent phase | |
| duration2 <- duration2[strike.id] | |
| ddif3 <- dmax | |
| ddif3$time <- abs(ddif3$time) | |
| ddif3$strike.num <- dmax$strike.num | |
| silent <- | |
| ddif3$time # silent phase = gap phase + artifact phase | |
| silent <- silent[strike.id] | |
| spike.num.over.duration <- | |
| spike.num / duration # discharge rate in co-activation | |
| spike.num.over.duration2 <- | |
| spike.num / duration2 # discharge rate including silent phase | |
| ## # spikes of initial 100 msec | |
| spike.init.time <- dmin$time | |
| spike.num.init100 <- c() | |
| for (k in 1:length(spike.init.time)) { | |
| spike.df2 <- spike.df[spike.df$time >= spike.init.time[k] & | |
| spike.df$time <= spike.init.time[k] + 0.1,] | |
| spike.df2 <- spike.df2[spike.df2$strike.num == k,] | |
| spike.num.init100 <- | |
| c(spike.num.init100, as.numeric(table(spike.df2$strike.num))) | |
| } | |
| spike.num.init100 <- spike.num.init100[strike.id] | |
| ## # spikes of last 100 msec in co-activation phase | |
| spike.last.time <- dmax$time | |
| spike.num.last100 <- c() | |
| for (k in 1:length(spike.last.time)) { | |
| spike.df3 <- spike.df[spike.df$time >= spike.last.time[k] - 0.1 & | |
| spike.df$time <= spike.last.time[k],] | |
| spike.df3 <- spike.df3[spike.df3$strike.num == k,] | |
| spike.num.last100 <- | |
| c(spike.num.last100, as.numeric(table(spike.df3$strike.num))) | |
| } | |
| spike.num.last100 <- spike.num.last100[strike.id] | |
| ## # spikes of last 100 msec including gap + artifact phase = silent phase | |
| spike.df4 <- | |
| spike.df[spike.df$time >= -0.1 & spike.df$time <= 0,] | |
| spike.num.100.with.silent <- | |
| as.numeric(table(spike.df4$strike.num)) | |
| spike.num.100.with.silent <- | |
| spike.num.100.with.silent[strike.id] | |
| cat("id:", list(...)[[i]]$id, "strike.num:", length(strike.id), "\n") | |
| df <- rbind( | |
| df, | |
| data.frame( | |
| id = list(...)[[i]]$id, | |
| ## kinematic variables | |
| ## meral-V kinematics | |
| strike.num = kinem.df$strike.num, | |
| input.rotation = kinem.df$input.rotation, # degrees | |
| input.rotation2 = kinem.df$input.rotation2, # degrees | |
| rotational.vel.m.in = kinem.df$rotational.vel.m.in / 180 * pi, # rad/s | |
| max.acc.m.in = kinem.df$max.acc.m.in / 180 * pi, # rad/s^2 | |
| velocity.meralV = kinem.df$velocity.meralV / 180 * pi, # rad/s | |
| max.acc.meralV = kinem.df$max.acc.meralV / 180 * pi, # rad/s^2 | |
| time.at.max.acc.meralV = kinem.df$time.at.max.acc.meralV * 10^3, # ms | |
| time.diff.max.acc = kinem.df$time.at.max.acc - kinem.df$time.at.max.acc.meralV, # ms | |
| ## strike duration | |
| time.at.limit = kinem.df$time.at.limit * 10 ^ 3, # ms | |
| time.at.max.linear.vel = kinem.df$time.at.max.linear.vel, | |
| time.at.max.acc = kinem.df$time.at.max.acc * 10 ^ 3, # ms | |
| ## strike linear kinematics | |
| max.linear.vel = kinem.df$max.linear.vel, | |
| linear.vel.at.limit = kinem.df$linear.vel.at.limit, | |
| max.linear.acc = kinem.df$max.linear.acc, | |
| linear.acc.at.limit = kinem.df$linear.acc.at.limit, | |
| ## strike angular kinematics | |
| max.vel = kinem.df$max.vel / 180 * pi, # rad/s | |
| vel.at.limit = kinem.df$vel.at.limit / 180 * pi, # rad/s | |
| max.acc = kinem.df$max.acc / 180 * pi, | |
| acc.at.limit = kinem.df$acc.at.limit / 180 * pi, | |
| ## EMG variables | |
| spike.num = spike.num, | |
| spike.num.init100 = spike.num.init100, | |
| spike.num.last100 = spike.num.last100, | |
| spike.num.100.with.silent = spike.num.100.with.silent, | |
| duration = duration, | |
| duration2 = duration2, | |
| silent = silent, # silent phase (gap + artifact) | |
| spike.num.over.duration = spike.num.over.duration, | |
| spike.num.over.duration2 = spike.num.over.duration2 # including silent phase | |
| ) | |
| ) | |
| } | |
| return(df) | |
| } | |
| aggregate.stats <- | |
| function(stat.df = stat.df, variable = "max.linear.vel", my.times = 10 ^ (-3), my.round = 2) { | |
| ## to generate descriptive statistics of stat.df | |
| dd <- stat.df[c("id", variable)] | |
| dd <- dd[complete.cases(dd),] | |
| mean.df <- aggregate(dd[-1], dd[1], mean, na.rm = T) | |
| mean.df[c(-1)] <- round(mean.df[c(-1)] * my.times, my.round) | |
| sd.df <- aggregate(dd[-1], dd[1], sd, na.rm = T) | |
| sd.df[c(-1)] <- round(sd.df[c(-1)] * my.times, my.round) | |
| range.df <- aggregate(dd[-1], dd[1], range, na.rm = T) | |
| range.df[,2] <- round(range.df[,2] * my.times, my.round) | |
| for (i in unique(mean.df$id)) { | |
| cat(i, ", ") | |
| cat(mean.df[mean.df$id == i, 2], " x ", paste(1 / my.times), "+-", | |
| sd.df[sd.df$id == i, 2], " x ", paste(1 / my.times), "\n") | |
| cat(i, ", ") | |
| cat( | |
| range.df[mean.df$id == i, 2][,1], " x ", paste(1 / my.times), " to ", | |
| range.df[mean.df$id == i, 2][,2], " x ", paste(1 / my.times), "\n" | |
| ) | |
| } | |
| } | |
| model.selection <- function(df = stat.df){ | |
| ## to perform model selection and output the AICs and coefficients | |
| df <- df[!is.na(df$strike.num),] | |
| res.var <- c("input.rotation", "vel.at.limit") | |
| exp.var <- c("input.rotation", "spike.num", "spike.num.init100", "spike.num.last100", | |
| "spike.num.100.with.silent", "duration", "silent", | |
| "spike.num.over.duration", "spike.num.over.duration2") | |
| exp.var.for.input <- exp.var[-1] | |
| df.input <- data.frame() | |
| for (i in 1:length(exp.var.for.input)) { | |
| my.formulas.input <- formula(paste(res.var[1], "~", exp.var.for.input[i])) | |
| my.formulas.input.random <- formula(paste("~1+", exp.var.for.input[i], "| id")) | |
| LMM.REML.input <- lme( my.formulas.input, | |
| random = my.formulas.input.random, | |
| data = df, control = lmeControl(opt = 'optim')) | |
| LMM.ML.input <- lme( my.formulas.input, | |
| random = my.formulas.input.random, | |
| data = df, control = lmeControl(opt = 'optim'), method = "ML") | |
| df.input <- rbind(df.input, data.frame(explanatory.var = paste(my.formulas.input)[3], | |
| AIC.REML = AIC(LMM.REML.input), | |
| AIC.ML = AIC(LMM.ML.input), | |
| intercept = fixed.effects(LMM.REML.input)[1], | |
| intercept.se = summary(LMM.REML.input)$tTable[1,2], | |
| slope = fixed.effects(LMM.REML.input)[2], | |
| slope.se = summary(LMM.REML.input)$tTable[2,2], | |
| pval = summary(LMM.REML.input)$tTable[2,5])) | |
| } | |
| df.output <- data.frame() | |
| for (i in 1:length(exp.var)) { | |
| my.formulas.output <- formula(paste(res.var[2], "~", exp.var[i])) | |
| my.formulas.output.random <- formula(paste("~1+", exp.var[i], "| id")) | |
| LMM.REML.output <- lme( my.formulas.output, | |
| random = my.formulas.output.random, | |
| data = df, control = lmeControl(opt = 'optim')) | |
| LMM.ML.output <- lme( my.formulas.output, | |
| random = my.formulas.output.random, | |
| data = df, control = lmeControl(opt = 'optim'), method="ML") | |
| df.output <- rbind(df.output, data.frame(explanatory.var = paste(my.formulas.output)[3], | |
| AIC.REML = AIC(LMM.REML.output), | |
| AIC.ML = AIC(LMM.ML.output), | |
| intercept = fixed.effects(LMM.REML.output)[1], | |
| intercept.se = summary(LMM.REML.output)$tTable[1,2], | |
| slope = fixed.effects(LMM.REML.output)[2], | |
| slope.se = summary(LMM.REML.output)$tTable[2,2], | |
| pval = summary(LMM.REML.output)$tTable[2,5])) | |
| } | |
| ## save the results of no random effects models | |
| LM.REML.input <- gls(input.rotation ~ spike.num, data = df, method = "REML") | |
| LM.REML.output <- gls(vel.at.limit ~ spike.num, data = df, method = "REML") | |
| ## combine the results no random effects models and other results | |
| df.input <- rbind(df.input, data.frame(explanatory.var = "no.random.effects", | |
| AIC.REML = AIC(LM.REML.input), | |
| AIC.ML = NA, | |
| intercept = coef(LM.REML.input)[1], | |
| intercept.se = summary(LM.REML.input)$tTable[1,2], | |
| slope = coef(LM.REML.input)[2], | |
| slope.se = summary(LM.REML.input)$tTable[2,2], | |
| pval = summary(LM.REML.input)$tTable[2,4])) | |
| df.output <- rbind(df.output, data.frame(explanatory.var = "no.random.effects", | |
| AIC.REML = AIC(LM.REML.output), | |
| AIC.ML = NA, | |
| intercept = coef(LM.REML.output)[1], | |
| intercept.se = summary(LM.REML.output)$tTable[1,2], | |
| slope = coef(LM.REML.output)[2], | |
| slope.se = summary(LM.REML.output)$tTable[2,2], | |
| pval = summary(LM.REML.output)$tTable[2,4])) | |
| ## null models | |
| Null.ML.input <- lme(input.rotation ~ 1, random = ~1|id, data = df, | |
| control = lmeControl(opt = 'optim'), method = "ML") | |
| Null.REML.input <- lme(input.rotation ~ 1, random = ~1|id, data = df, | |
| control = lmeControl(opt = 'optim'), method = "REML") | |
| Null.ML.output <- lme(vel.at.limit ~ 1, random = ~1|id, data = df, | |
| control = lmeControl(opt = 'optim'), method = "ML") | |
| Null.REML.output <- lme(vel.at.limit ~ 1, random = ~1|id, data = df, | |
| control = lmeControl(opt = 'optim'), method = "REML") | |
| df.input <- rbind(df.input, data.frame(explanatory.var = "null", | |
| AIC.REML = AIC(Null.REML.input), | |
| AIC.ML = AIC(Null.ML.input), | |
| intercept = fixed.effects(Null.REML.input)[1], | |
| intercept.se = summary(Null.REML.input)$tTable[1,2], | |
| slope = 0, | |
| slope.se = NA, | |
| pval = NA)) | |
| df.output <- rbind(df.output, data.frame(explanatory.var = "null", | |
| AIC.REML = AIC(Null.REML.output), | |
| AIC.ML = AIC(Null.ML.output), | |
| intercept = fixed.effects(Null.REML.output)[1], | |
| intercept.se = summary(Null.REML.output)$tTable[1,2], | |
| slope = 0, | |
| slope.se = NA, | |
| pval = NA)) | |
| ## print the tables | |
| print(df.input) | |
| print(df.output) | |
| } | |
| fit.input.rotation.vs.number.of.spike <- | |
| function(df = stat.df) { | |
| ## to explain the relationship of input spring rotation and the number of spikes (Fig.10A) | |
| ddf <- df[!is.na(df$strike.num),] | |
| ## df must be generated using generate.data.for.model.selection | |
| Mlme1 <- lme(input.rotation ~ 1, random = ~ 1 | id, data = ddf) | |
| Mlme2 <- lme( | |
| input.rotation ~ spike.num, random = ~ 1 + spike.num | id, | |
| control = lmeControl(opt = 'optim'), data = ddf | |
| ) | |
| ddf$M1F0 <- fitted(Mlme1, level = 0) | |
| ddf$M1F1 <- fitted(Mlme1, level = 1) | |
| ddf$M2F0 <- fitted(Mlme2, level = 0) | |
| ddf$M2F1 <- fitted(Mlme2, level = 1) | |
| graphics.off() | |
| my.plot1 <- ggplot(ddf) + | |
| geom_point(aes( | |
| spike.num, input.rotation, shape = id, colour = id | |
| )) + | |
| geom_line( | |
| aes(spike.num, M1F0), lwd = 2, col = "grey", alpha = 0.2, lty = 2 | |
| ) + | |
| geom_line(aes(spike.num, M2F0), lwd = 2, alpha = 0.5, col = "grey") + | |
| geom_line(aes(spike.num, M1F1, group = id, colour = id), alpha = 0.2, lty = 2) + | |
| geom_line(aes(spike.num, M2F1, group = id, colour = id), alpha = 0.8) + | |
| scale_color_manual(values = c("blue", "dark orange", | |
| "red", "purple", "dark green")) + | |
| scale_shape_manual(values = c(1,2,3,4,5)) + | |
| xlab("# spikes") + | |
| ylab("net meral-V rotation (degree)") + | |
| my.theme() | |
| print(my.plot1) | |
| cat("* The summary of the constant model:\n") | |
| print(summary(Mlme1)) | |
| cat("\n") | |
| cat("* The summary of the correlative model:\n") | |
| print(summary(Mlme2)) | |
| } | |
| fit.strike.velocity.vs.number.of.spike <- function(df = stat.df) { | |
| ## to explain the relationship of strike velocity and the number of spikes (Fig.10B) | |
| df <- df[!is.na(df$strike.num),] | |
| ## df must be generated using generate.data.for.model.selection | |
| Mlme1 <- lme(vel.at.limit ~ 1, random = ~ 1 | id, data = df) | |
| Mlme2 <- lme( | |
| vel.at.limit ~ spike.num, random = ~ 1 + spike.num | id, | |
| control = lmeControl(opt = 'optim'), data = df | |
| ) | |
| df$M1F0 <- fitted(Mlme1, level = 0) | |
| df$M1F1 <- fitted(Mlme1, level = 1) | |
| df$M2F0 <- fitted(Mlme2, level = 0) | |
| df$M2F1 <- fitted(Mlme2, level = 1) | |
| graphics.off() | |
| my.plot1 <- ggplot(df) + | |
| geom_point(aes(spike.num, vel.at.limit, shape = id, colour = id)) + | |
| geom_line( | |
| aes(spike.num, M1F0), lwd = 2, col = "grey", alpha = 0.2, lty = 2) + | |
| geom_line(aes(spike.num, M2F0), | |
| lwd = 2, alpha = 0.5, col = "grey") + | |
| geom_line(aes(spike.num, | |
| M1F1, group = id, colour = id), alpha = 0.2, lty = 2) + | |
| geom_line(aes(spike.num, | |
| M2F1, group = id, colour = id), alpha = 0.8) + | |
| scale_color_manual( | |
| values = c("blue", "dark orange", "red", "purple", "dark green")) + | |
| scale_shape_manual(values = c(1,2,3,4,5)) + | |
| xlab("# spikes") + | |
| ylab("strike velocity (rad/sec)") + | |
| my.theme() | |
| print(my.plot1) | |
| cat("* The summary of the constant model:\n") | |
| print(summary(Mlme1)) | |
| cat("\n") | |
| cat("* The summary of the correlative model:\n") | |
| print(summary(Mlme2)) | |
| } | |
| plot.time.diff.acc <- function(d){ | |
| ## to plot time difference (Fig.8) | |
| ggplot(d) + geom_point(aes(strike.num, time.diff.max.acc*1000, size = vel.at.limit), pch = 1) + | |
| geom_hline(aes(yintercept = 0), linetype = 2) + | |
| xlab("Strike number") + ylab("Time interval between max. acceleration \n of propodus and meral-V (ms)") + | |
| facet_grid(~id) + my.theme() | |
| } | |
| plot.emg.contraction <- | |
| function(d = nb30, emg.ch = 1, strike.num = 1) { | |
| ## to plot input meral-V rotation and EMG of the lateral extensor (Fig.7A) | |
| emg <- d$phys$aligned[[strike.num]]$phys | |
| rot <- -d$contraction$input.rotation[[strike.num]] | |
| fps <- d$fps.contraction | |
| x <- -seq(0, 1 / fps * (length(rot) - 1), 1 / fps) | |
| x <- x[length(x):1] | |
| emg <- emg[, c(1, emg.ch + 1)] | |
| colnames(emg) <- c("time", "value") | |
| emg <- subset(emg, time >= -0.9 & time <= 0) | |
| start.time <- min(x) | |
| emg$value <- emg$value * 10 | |
| rot.df <- data.frame(time = x, value = rot) | |
| rot.df <- | |
| data.frame(approx(rot.df$time, rot.df$value, xout = emg$time)) | |
| colnames(rot.df) <- c("time", "value") | |
| rot.df <- cbind(rot.df, category = "rotation") | |
| emg <- cbind(emg, category = "emg") | |
| df <- rbind(rot.df, emg) | |
| p <- qplot(time, value, data = df, geom = "path") + | |
| facet_grid(category ~ .) + my.theme() | |
| return(p) | |
| } | |
| fit.3.kinematics <- function(d = nb30, my.strikes = c(1,2,3)) { | |
| ## to fit the 10th order polinomial and GAM fittings of 3 strikes (Fig.7B) | |
| d1 <- fit.kinematics(d, my.strikes[1]) | |
| d2 <- fit.kinematics(d, my.strikes[2]) | |
| d3 <- fit.kinematics(d, my.strikes[3]) | |
| mrd <- rbind(d1$mrd, d2$mrd, d3$mrd) | |
| prd <- rbind(d1$prd, d2$prd, d3$prd) | |
| fit.df <- rbind(d1$fit.df, d2$fit.df, d3$fit.df) | |
| graphics.off() | |
| print( ggplot(mrd) + | |
| geom_point(aes(time * 10 ^ 3,-rotation), alpha = 0.8, pch = 1) + | |
| geom_line( aes(time * 10 ^ 3,-pred.poly.m), col = 2, alpha = 0.5, data = fit.df ) + | |
| geom_line( aes(time * 10 ^ 3,-pred.gam.m), alpha = 0.5, data = fit.df ) + | |
| geom_point( aes(time * 10 ^ 3, rotation), alpha = 0.8, pch = 1, data = prd ) + | |
| geom_line( aes(time * 10 ^ 3, pred.poly.p), col = 2, alpha = 0.5, data = fit.df ) + | |
| geom_line( aes(time * 10 ^ 3, pred.gam.p, group = strike.num), alpha = 0.5, data = fit.df ) + | |
| geom_hline(aes(yintercept = 20), lty = 2) + | |
| geom_hline(aes(yintercept = 0)) + | |
| xlab("time (msec)") + ylab("rotation (degree)") + | |
| facet_grid(strike.num ~ .) + | |
| my.theme() | |
| ) | |
| } | |
| fit.kinematics <- function(d = nb30, my.strike.num = 1) { | |
| ## to fit the 10th order polinomial and GAM fittings of one strikes (Fig.7B) | |
| mr <- d$strike$output.meralV.rotation | |
| pr <- d$strike$output.propodus.rotation | |
| fps <- d$fps.strike | |
| mrd <- data.frame() | |
| prd <- data.frame() | |
| for (i in 1:length(mr)) { | |
| if (length(mr[[i]]) > 0) { | |
| mrd <- rbind(mrd, data.frame( | |
| strike.num = i, | |
| rotation = mr[[i]], | |
| time = seq(0, 1 / fps * (length(mr[[i]]) - 1), 1 / fps) | |
| )) | |
| } else { | |
| mrd <- rbind(mrd, data.frame( | |
| strike.num = i, rotation = NA, time = NA | |
| )) | |
| } | |
| if (length(pr[[i]]) > 0) { | |
| prd <- rbind(prd, data.frame( | |
| strike.num = i, | |
| rotation = pr[[i]], | |
| time = seq(0, 1 / fps * (length(pr[[i]]) - 1), 1 / fps) | |
| )) | |
| } else { | |
| prd <- rbind(prd, data.frame( | |
| strike.num = i, rotation = NA, time = NA | |
| )) | |
| } | |
| } | |
| mrd <- subset(mrd, strike.num == my.strike.num) | |
| prd <- subset(prd, strike.num == my.strike.num) | |
| max.len <- min(dim(mrd)[1], dim(prd)[1]) | |
| mrd <- mrd[1:max.len,] | |
| prd <- prd[1:max.len,] | |
| ## gam fitting | |
| fit.gam.m <- gam(rotation ~ s(time, k = 5), data = mrd) | |
| time.grid <- seq(min(mrd$time), max(mrd$time), 1 / fps / 10) | |
| pred.gam.m <- | |
| predict(fit.gam.m, newdata = data.frame(time = time.grid), type = "response") | |
| fit.gam.p <- gam(rotation ~ s(time, k = 5), data = prd) | |
| time.grid <- seq(min(prd$time), max(prd$time), 1 / fps / 10) | |
| pred.gam.p <- | |
| predict(fit.gam.p, newdata = data.frame(time = time.grid), type = "response") | |
| ## polynomial fitting | |
| mrd.poly <- mrd | |
| for (i in 1:10) { | |
| mrd.poly <- rbind(mrd[1,], mrd.poly) | |
| } | |
| for (i in 1:10) { | |
| mrd.poly <- rbind(mrd.poly, mrd[dim(mrd)[1],]) | |
| } | |
| mrd.poly$time <- 0:(dim(mrd.poly)[1] - 1) | |
| prd.poly <- prd | |
| for (i in 1:10) { | |
| prd.poly <- rbind(prd[1,], prd.poly) | |
| } | |
| for (i in 1:10) { | |
| prd.poly <- rbind(prd.poly, prd[dim(prd)[1],]) | |
| } | |
| prd.poly$time <- 0:(dim(prd.poly)[1] - 1) | |
| fit.poly.m <- lm(rotation ~ poly(time, 10), data = mrd.poly) | |
| time.grid <- seq(min(mrd.poly$time), max(mrd.poly$time), 1 / 10) | |
| pred.poly.m <- | |
| predict(fit.poly.m, newdata = data.frame(time = time.grid), type = "response") | |
| for (i in 1:100) { | |
| pred.poly.m <- pred.poly.m[-1] | |
| } | |
| for (i in 1:100) { | |
| pred.poly.m <- pred.poly.m[-length(pred.poly.m)] | |
| } | |
| fit.poly.p <- lm(rotation ~ poly(time, 10), data = prd.poly) | |
| time.grid <- seq(min(prd.poly$time), max(prd.poly$time), 1 / 10) | |
| pred.poly.p <- | |
| predict(fit.poly.p, newdata = data.frame(time = time.grid), type = "response") | |
| for (i in 1:100) { | |
| pred.poly.p <- pred.poly.p[-1] | |
| time.grid <- time.grid[-1] | |
| } | |
| for (i in 1:100) { | |
| pred.poly.p <- pred.poly.p[-length(pred.poly.p)] | |
| time.grid <- time.grid[-length(time.grid)] | |
| } | |
| time.grid <- time.grid - time.grid[1] | |
| fit.df <- data.frame( | |
| strike.num = my.strike.num, | |
| time = time.grid / fps, pred.poly.m = pred.poly.m, pred.gam.m = | |
| pred.gam.m, | |
| pred.poly.p = pred.poly.p, pred.gam.p = pred.gam.p | |
| ) | |
| print( | |
| ggplot(mrd) + geom_point(aes(time * 10 ^ 3,-rotation), pch = 1) + | |
| geom_line( | |
| aes(time * 10 ^ 3,-pred.poly.m), col = 2, data = fit.df | |
| ) + | |
| geom_line(aes(time * 10 ^ 3,-pred.gam.m), data = fit.df) + | |
| geom_point(aes(time * 10 ^ 3, rotation), pch = 1, data = prd) + | |
| geom_line( | |
| aes(time * 10 ^ 3, pred.poly.p), col = 2, data = fit.df | |
| ) + | |
| geom_line(aes(time * 10 ^ 3, pred.gam.p), data = fit.df) + | |
| xlab("time (msec)") + ylab("rotation (degree)") + | |
| my.theme() | |
| ) | |
| return(list( | |
| strike.num = my.strike.num, mrd = mrd, prd = prd, fit.df = fit.df | |
| )) | |
| } | |
| my.theme <- function() { | |
| ## to custamize the ggplot for publication | |
| ## usage: q + my.theme() | |
| theme( | |
| panel.grid.major = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| panel.background = element_blank() | |
| ) + | |
| theme( | |
| axis.line = element_line(), | |
| axis.text = element_text(colour = "black"), | |
| axis.ticks = element_line(colour = "black") | |
| ) | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment