A method of compensating for the excitability of spinal motoneurones when estimating the magnitude of potentials evoked in skeletal muscles.

1. #!/usr/bin/Rscript --slave
3. # annotated_QR_example_code.R
4. # 5th May 2026
5. # Modified: 5th May 2026
7. # Author: Richard G. Carson ([email hidden])
9. #################
10. # Read the data file
11. data <- read.csv("example_data.csv", header=TRUE, stringsAsFactors=FALSE)
13. # Ensure data are numeric
14. data["mep_amplitude"] <- as.numeric(unlist(data["mep_amplitude"]))
15. data["rms_emg"] <- as.numeric(unlist(data["rms_emg"]))
17. ###################################
19. # Whereas least-squares regression (LSR) estimates the conditional mean as a linear
20. # combination of the predictors, quantile regression (QR) estimates a conditional quantile
21. # function (such as the 0.5 quantile, i.e., the median) as a linear combination of the
22. # predictors.
23. # Unlike least squares regression, QR does not entail an assumption that the outcome
24. # variable has constant variance or that samples conform to a particular parametric
25. # distribution. It is also robust to the presence of outliers. It was implemented here
26. # by means of the rq() function from the quantreg package (Koenker, 2022).
27. # The “tau” parameter was specified as 0.5, indicating that the conditional quantile
28. # for the median was to be used.
30. ###################################
31. library(quantreg)
33. # Estimate the quantile regression model for the median (tau = 0.5)
34. qr_model <- rq(mep_amplitude ~ rms_emg, data = data, tau = 0.5)
36. # Generate residuals
37. data$qr_residuals <- qr_model$residuals
39. ###################################
40. # A positive value of a residual indicates that an estimate of MEP amplitude is larger
41. # than that which would be predicted by r.m.s. EMG alone – to an extent corresponding
42. # to its magnitude. A negative residual value indicates an estimate of MEP amplitude
43. # that is smaller than would be predicted for the r.m.s. EMG recorded in the period
44. # preceding the stimulation.
45. # The use of the unadjusted residuals as the basis for further analysis is however
46. # limited by the fact that for a given sample the sum of the residuals necessarily tends
47. # to zero. It is therefore desirable that the residuals are expressed relative
48. # to an appropriate reference value.
49. # The parameters of the linear regression model used to generate each set of residuals
50. # affords a means of addressing this requirement.
51. ###################################
53. ###################################
54. # Extract key information from the fitted model
55. # y = mx + c
56. # y - c = mx
57. # (y-c)/m = x
58. # m = lm.mod_Ch5$coeff[["rms_emg"]]
59. # c = lm.mod_Ch5$coeff[["(Intercept)"]]
60. qr.mod_slope <- qr_model$coeff[["rms_emg"]]
61. qr.mod_intercept <- qr_model$coeff[["(Intercept)"]]
63. ###################################
65. ###################################
66. # In an idealised scenario in which a linear relationship between r.m.s. EMG and MEP
67. # amplitude accounts for 100% of the variance present in the data, the appropriate
68. # reference is obtained as the value of the intercept (i.e., with the ordinate (y) axis).
69. # This corresponds to the predicted value of the MEP amplitude when the r.m.s. EMG is
70. # equal to zero. In this idealised case, the adjusted MEP amplitudes for the sample are
71. # simply obtained by adding the value of the intercept to each of the residuals.
72. # If there is no relationship between r.m.s. EMG and MEP amplitude, the best estimate
73. # for all values of r.m.s. EMG is the central tendency of the original MEP amplitude
74. # values.
75. # That is. if the slope of the regression line is equal to zero, the reference value is
76. # provided by the central tendency of the actual values (which necessarily equates to the
77. # central tendency of the predicted values and, in this specific case only, to the value
78. # of the intercept).
79. # Thus, the extent to which the appropriate reference value differs from the central
80. # tendency of the actual values is contingent on the slope of the regression line.
82. # With respect to real data however, a linear model will not typically account for
83. # all variation inherent to the sample. The use of a linear regression model
84. # dictates that there is no uncertainty about the estimated slope at the central tendency
85. # of the “predictor” variable (the r.m.s. EMG in the present instance). Even though
86. # prediction will also be most precise at the estimate of central tendency, there will
87. # nonetheless be some uncertainty concerning the true vertical position of the regression
88. # line at this point.
89. # The degree of uncertainty increases with separation from the centre of the distribution
90. # of the predictor variable. That is, unless the value of the intercept coincides with
91. # the value of the estimate of central tendency, the standard error of the ordinate at
92. # the intercept will always be larger than the standard error of the ordinate at the
93. # central tendency of the “predictor” variable.
95. # Since the standard error of the estimate is a measure of the accuracy of prediction,
96. # we wish to accord to the intercept (in its contribution to the generation of a
97. # reference value) a weighting that varies inversely with its standard error.
98. # Specifically, as the standard error of the estimate at the intercept tends to infinity,
99. # the contribution of the intercept to the generation of a reference value should tend
100. # towards zero. As noted above, as the slope of the regression line approaches zero,
101. # the reference value derived from the intercept will converge upon that of the central
102. # tendency of the original values. Accordingly, therefore, the relative weighting of the
103. # intercept should also be such that, with increases in its standard error, the reference
104. # value converges on the central tendency. In summary, the extent to which the reference
105. # value used to adjust the residuals differs from the central tendency of the
106. # predicted/actual values, will depend on the slope of the regression line and the
107. # uncertainty with which its intercept can be estimated.
108. ###################################
110. ###################################
111. # Determine the value of x (i.e., the r.m.s. EMG) that corresponds to
112. # the median of the predicted values.
113. qr_predicted_values <- fitted(qr_model)
114. qr_predicted_median <- median(qr_predicted_values)
115. qr_rms_emg_at_median <- (qr_predicted_median - qr_model$coeff[["(Intercept)"]] )/qr_model$coeff[["rms_emg"]]
117. # Determine the "quasi" standard error of the fit when the y value (MEP amplitude)
118. # is predicted by means of the qr model when the x value (r.m.s. EMG)
119. # is that at which the median of the predicted values is obtained.
120. qr_predict_at_median <- predict(qr_model, newdata=data.frame(rms_emg = qr_rms_emg_at_median), interval = "confidence")
121. qr_predict_at_median_ci <- qr_predict_at_median[1,"higher"] - qr_predict_at_median[1,"lower"]
123. # As confidence intervals (i.e., 95%) for prediction are reported by the rq() function
124. # rather than standard errors, the values were first converted to the equivalent
125. # standard errors
127. qr_predict_at_median_SE <- qr_predict_at_median_ci / qnorm(0.975)
129. # Determine the "quasi" standard error of the fit when the y value (MEP amplitude)
130. # is predicted by means of the linear model when the x value (r.m.s. EMG)
131. # is equal to zero, i.e. the intercept.
132. qr_predict_at_intercept <- predict(qr_model, newdata=data.frame(rms_emg=0.0), interval = "confidence")
133. qr_predict_at_intercept_ci <- qr_predict_at_intercept[1,"higher"] - qr_predict_at_intercept[1,"lower"]
135. # As confidence intervals (i.e., 95%) for prediction are reported by the rq() function
136. # rather than standard errors, the values were first converted to the equivalent
137. # standard errors
139. qr_predict_at_intercept_SE <- qr_predict_at_intercept_ci / qnorm(0.975)
141. # Calculate the sum of the SE for the predicted value of y at the median and at the
142. # intercept.
143. qr_predict_at_median_AND_intercept_SE_sum <- qr_predict_at_median_SE + qr_predict_at_intercept_SE
145. # Express the SE for the predicted value of y at the median and at the intercept as a
146. # proportion of the sum.
147. qr_predict_at_median_proportion <- qr_predict_at_median_SE / qr_predict_at_median_AND_intercept_SE_sum
148. qr_predict_at_intercept_proportion <- qr_predict_at_intercept_SE / qr_predict_at_median_AND_intercept_SE_sum
150. # Calculate median weighting and weighted values
151. qr_median_weight <- (qr_predict_at_intercept_SE/qr_predict_at_median_AND_intercept_SE_sum) - (qr_predict_at_median_SE/qr_predict_at_median_AND_intercept_SE_sum)
152. qr_median_weighted_value <- qr_median_weight * qr_predicted_median
154. # Calculate intercept weighting and weighted values
155. qr_intercept_weight <- (1 - (qr_predict_at_intercept_SE/qr_predict_at_median_AND_intercept_SE_sum)) + (qr_predict_at_median_SE/qr_predict_at_median_AND_intercept_SE_sum)
156. qr_intercept_weighted_value <- qr_intercept_weight * qr.mod_intercept
158. # Generate final estimate of the reference value
159. qr_weighted_reference_amplitude <- qr_median_weighted_value + qr_intercept_weighted_value
161. # Add the final estimate to the residuals generated using the quantile regression
162. data$qr_adjusted_mep_amplitude <- qr_model$residuals + qr_weighted_reference_amplitude
164. detach("package:quantreg")
166. ###################################
169. ###################################
170. # Write the data to an .xlsx format file
171. ###################################
173. # load the openxlsx package
174. library(openxlsx)
176. # In order to permit data to be written to an .xlsx file, create a workbook
177. wbOut<-createWorkbook()
179. # Add worksheets to the workbook object
180. addWorksheet(wbOut, "adjusted")
182. # Write the summarised data to the workbook object
183. writeData(wbOut, "adjusted", data, startCol = 1, startRow = 1, xy = NULL, colNames = TRUE, rowNames = FALSE, headerStyle = NULL, keepNA = TRUE)
185. # Save the summary data workbook - overwriting the original file (if it exists)
186. saveWorkbook(wbOut, "example_data_qr_adjusted.xlsx", overwrite = TRUE)
188. detach("package:openxlsx")
190. ###################################

annotated_QR_example_code.R, under CC-BY-4.0 · at the source