Skip to content

Commit 06d1efe

Browse files
upd spacing
1 parent d0fb02c commit 06d1efe

File tree

4 files changed

+685
-484
lines changed

4 files changed

+685
-484
lines changed

R/01-R_PythonSetup.Rmd

Lines changed: 7 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ This R Markdown file corresponds to data analysis for a **Developmental Cognitiv
88

99
This tutorial shows how to access Python in R Studio via the *reticulate* package (Ushey et al., 2020). While there are multiple ways in which a user may call Python from R (e.g., **repl_python()** function), our discussion focuses solely on how to access Python in an R Markdown file. Please see the [*reticulate* website](https://rstudio.github.io/reticulate) for details on alternative approaches. We also recommend reviewing the [reticulate cheat sheet](https://ugoproto.github.io/ugo_r_doc/pdf/reticulate.pdf) provided by R Studio.
1010

11-
## **1.** Install packages
11+
1212
### Load reticulate
1313
```{r}
1414
#install.packages("reticulate")
@@ -27,15 +27,14 @@ Several Python (and conda) versions may already exist on your computer. It is th
2727
Sys.which("python")
2828
```
2929

30-
## **2.** Identify and set python
3130
### Return all detected python versions
3231
```{r}
3332
py_discover_config()
3433
```
3534

3635
### List all available conda environments
3736
```{r}
38-
conda_list(conda= "auto")
37+
conda_list(conda = "auto")
3938
```
4039

4140
### Set python version
@@ -50,14 +49,12 @@ usethis::edit_r_profile()
5049
use_virtualenv("r-reticulate")
5150
```
5251

53-
54-
## **3.** Load python packages
5552
### Install packages
5653
Uncomment this code when you are ready to install python packages.
5754
```{r}
58-
#py_install("pandas", envname= "r-reticulate")
59-
#py_install("numpy", envname= "r-reticulate")
60-
#py_install("fooof", envname= "r-reticulate")
61-
#py_install("matplotlib", envname= "r-reticulate")
62-
#py_install("pathlib", envname= "r-reticulate")
55+
#py_install("pandas", envname = "r-reticulate")
56+
#py_install("numpy", envname = "r-reticulate")
57+
#py_install("fooof", envname = "r-reticulate")
58+
#py_install("matplotlib", envname = "r-reticulate")
59+
#py_install("pathlib", envname = "r-reticulate")
6360
```

R/02-R_IndividualPSD.Rmd

Lines changed: 49 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -11,16 +11,16 @@ Results are saved out to a **Output** folder for further consideration. Please n
1111

1212
### Load R packages
1313
These packages are required for this tutorial. Install packages using the *install.packages()* function if needed.
14-
```{r, message= FALSE}
15-
library(tidyverse) # Access a collection of packages for data management (e.g., dplyr, ggplot2, readr)
14+
```{r, message = FALSE}
15+
library(tidyverse) # Access a collection of packages for data management
1616
library(reticulate) # Interface Python and R Studio
1717
library(magick) # Load and adjust .PNG files
1818
library(gridExtra) # Arrange multiple plots
1919
```
2020

2121
### Load Python libraries
2222
Load the Python modules we need for this tutorial.
23-
```{python, message=FALSE}
23+
```{python, message = FALSE}
2424
# Import some useful standard library modules
2525
import os
2626
from pathlib import Path
@@ -42,7 +42,7 @@ from fooof.analysis.error import compute_pointwise_error_fm
4242
```
4343

4444
### Check specparam version
45-
```{python, message=FALSE}
45+
```{python, message = FALSE}
4646
import fooof
4747
print(fooof.__version__)
4848
```
@@ -51,42 +51,42 @@ print(fooof.__version__)
5151
Load CSV files, including:
5252
- `freqs.csv`, which contains a vector of frequencies
5353
- `eopPSDs.csv`, which contains the power values for an individual power spectrum
54-
```{python, message=FALSE}
54+
```{python, message = FALSE}
5555
# Load csv files containing frequency and power values
5656
freqs = np.ravel(pd.read_csv("../Data/freqs.csv"))
5757
spectrum = np.ravel(pd.read_csv("../Data/indv.csv"))[1:101]
5858
```
5959

60-
```{python, message=FALSE}
60+
```{python, message = FALSE}
6161
# Check shape of loaded data
6262
print(freqs.shape)
6363
print(spectrum.shape)
6464
```
6565

6666
### Parameterize a power spectrum
6767
Now we can parameterize our power spectrum!
68-
```{python, message=FALSE}
68+
```{python, message = FALSE}
6969
# Initialize a model object for spectral parameterization, with some settings
70-
fm = FOOOF(peak_width_limits = [2, 8], max_n_peaks = 6, min_peak_height = 0.05, verbose = False)
70+
fm = FOOOF(peak_width_limits = [1, 8], max_n_peaks = 6, min_peak_height = 0.10, verbose = False)
7171
```
7272

73-
```{python, message=FALSE, fig.height= 6, fig.width= 9}
73+
```{python, message = FALSE, fig.height = 6, fig.width = 9}
7474
# Fit individual PSD over 3-40 Hz range
7575
fm.report(freqs, spectrum, [3, 40])
7676
plt.show()
7777
```
7878

79-
```{python, message=FALSE}
79+
```{python, message = FALSE}
8080
# Save out a copy of the model fit report
81-
fm.save_report('INDV_demo', file_path= '../Output')
81+
fm.save_report('INDV_demo', file_path = '../Output')
8282
8383
# The following line can also be used to save out the model plot
8484
#fm.plot(save_fig = True, file_name = 'INDV_demo', file_path = '../Output')
8585
```
8686

8787
### Access model fit information
8888
All of the model fit information is saved to the model object. Note that all attributes learned in the model fit process have a trailing underscore.
89-
```{python, message=FALSE}
89+
```{python, message = FALSE}
9090
# Access the model fit parameters from the model object
9191
print('Aperiodic parameters: \n', fm.aperiodic_params_, '\n')
9292
print('Peak parameters: \n', fm.peak_params_, '\n')
@@ -99,46 +99,46 @@ print('Number of fit peaks: \n', fm.n_peaks_)
9999

100100
### Extract periodic and aperiodic parameters
101101
Another way to access model fit parameters is to use the `get_params` method, which can be used to access model fit attributes.
102-
```{python, message=FALSE}
102+
```{python, message = FALSE}
103103
# Extract aperiodic and periodic parameter information
104104
aps = fm.get_params('aperiodic_params')
105105
peaks = fm.get_params('peak_params')
106106
```
107107

108-
```{python, message=FALSE}
108+
```{python, message = FALSE}
109109
# Extract goodness of fit information
110110
err = fm.get_params('error')
111111
r2s = fm.get_params('r_squared')
112112
```
113113

114-
```{python, message=FALSE}
114+
```{python, message = FALSE}
115115
# Extract specific parameters
116116
exp = fm.get_params('aperiodic_params', 'exponent')
117117
cfs = fm.get_params('peak_params', 'CF')
118118
```
119119

120-
```{python, message=FALSE}
120+
```{python, message = FALSE}
121121
# Print out a custom parameter report
122122
template = ("With an error level of {error:1.2f}, specparam fit an exponent of {exponent:1.2f} and peaks of {cfs:s} Hz.")
123123
print(template.format(error = err, exponent = exp, cfs = ' & '.join(map(str, [round(CF, 2) for CF in cfs]))))
124124
```
125125

126126
### Plot flattened power spectrum
127127
It may be useful to plot a flattened power spectrum, with the aperiodic fit removed.
128-
```{python, message=FALSE}
128+
```{python, message = FALSE}
129129
# Set whether to plot in log-log space
130130
plt_log = False
131131
```
132132

133-
```{python, message=FALSE}
133+
```{python, message = FALSE}
134134
# Do an initial aperiodic fit - a robust fit, that excludes outliers
135135
init_ap_fit = gen_aperiodic(fm.freqs, fm._robust_ap_fit(fm.freqs, fm.power_spectrum))
136136
137137
# Recompute the flattened spectrum using the initial aperiodic fit
138138
init_flat_spec = fm.power_spectrum - init_ap_fit
139139
```
140140

141-
```{python, message=FALSE}
141+
```{python, message = FALSE}
142142
# Plot the flattened the power spectrum
143143
plot_spectrum(fm.freqs, init_flat_spec, plt_log, label = 'Flattened spectrum', color = 'black')
144144
plt.show()
@@ -151,59 +151,70 @@ Transferring specparam model objects into R allows a researcher to work within a
151151
```{r}
152152
# Transfer periodic parameters to R data frame
153153
per <- as.data.frame(py$peaks) %>%
154-
rename(CF= 1, PW= 2, BW= 3) %>%
155-
mutate(peak_num= row_number()) %>% group_by(peak_num) %>%
156-
mutate(index= seq_along(CF)) %>%
157-
pivot_wider(id_col= index, names_from= peak_num, values_from= c(CF,PW,BW)) %>%
158-
select(index,CF_1,PW_1,BW_1,CF_2,PW_2,BW_2)
154+
rename(CF = 1, PW = 2, BW = 3) %>%
155+
mutate(peak_num = row_number()) %>%
156+
group_by(peak_num) %>%
157+
mutate(index = seq_along(CF)) %>%
158+
pivot_wider(id_col = index, names_from = peak_num, values_from = c(CF, PW, BW)) %>%
159+
select(index, CF_1, PW_1, BW_1, CF_2, PW_2, BW_2)
159160
```
160161

161162
```{r}
162163
# Transfer aperiodic parameters to R data frame
163-
aps <- as.data.frame(py$aps) %>% mutate(var= c("offset","exponent")) %>% spread(var,`py$aps`) %>% mutate(index= row_number()) %>% select(index,offset,exponent)
164+
aps <- as.data.frame(py$aps) %>%
165+
mutate(var = c("offset","exponent")) %>%
166+
spread(var, `py$aps`) %>% mutate(index = row_number()) %>%
167+
select(index, offset, exponent)
164168
```
165169

166170
```{r}
167171
# Transfer group fit information to R data frame
168-
r2s <- as.data.frame(py$r2s) %>% rename(r2s= 1) %>% mutate(index= row_number())
169-
err <- as.data.frame(py$err) %>% rename(err= 1) %>% mutate(index= row_number())
172+
r2s <- as.data.frame(py$r2s) %>%
173+
rename(r2s = 1) %>%
174+
mutate(index = row_number())
175+
176+
err <- as.data.frame(py$err) %>%
177+
rename(err = 1) %>%
178+
mutate(index = row_number())
170179
```
171180

172181
```{r}
173182
# Pull IDs
174-
IDs <- read.csv("../Data/indv.csv", header= TRUE) %>% select(ID) %>% mutate(index= row_number())
183+
IDs <- read.csv("../Data/indv.csv", header = TRUE) %>%
184+
select(ID) %>%
185+
mutate(index = row_number())
175186
176187
# Join data frames
177-
dat <- full_join(IDs, per, by= "index") %>%
178-
full_join(aps, by= "index") %>%
179-
full_join(r2s, by= "index") %>%
180-
full_join(err, by= "index") %>%
188+
dat <- full_join(IDs, per, by = "index") %>%
189+
full_join(aps, by = "index") %>%
190+
full_join(r2s, by = "index") %>%
191+
full_join(err, by = "index") %>%
181192
select(-index) %>% arrange(ID)
182193
```
183194

184195
### Save fit information
185196
```{r}
186197
# Save out data frame as csv file
187-
write.csv(dat, "../Output/INDV_demo.csv", row.names= FALSE)
198+
write.csv(dat, "../Output/INDV_demo.csv", row.names = FALSE)
188199
```
189200

190201
### Frequency-by-frequency error
191202
It can be useful to plot frequency-by-frequency error of the model fit, to identify where in frequency space the spectrum is (or is not) being fit well. When fitting individual spectrum, this can be accomplished using the `compute_pointwise_error_fm` function.
192203

193-
In this case, we can see that error fluctuates around 0.05, which is the same as the mean absolute error for the model (MAE). There are points in the spectrum where the model fit is somewhat poor, particularly < 4 Hz, ~6-9 Hz, and ~14 Hz. Further considerations may be necessary for this model fit.
194-
```{python, message=FALSE, fig.height= 6, fig.width= 9}
204+
In this case, we can see that error fluctuates around 0.06, which is the mean absolute error for the model (MAE). There are points in the spectrum where the model fit is somewhat poor, particularly < 4 Hz, ~6-9 Hz, and ~14 Hz. Further considerations may be necessary for this model fit.
205+
```{python, message = FALSE, fig.height = 6, fig.width = 9}
195206
# Plot frequency-by-frequency error
196207
compute_pointwise_error_fm(fm, plot_errors = True)
197208
plt.show()
198209
```
199210

200-
```{python, message=FALSE}
211+
```{python, message = FALSE}
201212
# Return the frequency-by-frequency errors
202-
errs_fm = compute_pointwise_error_fm(fm, plot_errors= False, return_errors = True)
213+
errs_fm = compute_pointwise_error_fm(fm, plot_errors = False, return_errors = True)
203214
errs_fm
204215
```
205216

206-
```{python, message=FALSE}
217+
```{python, message = FALSE}
207218
# Note that the average of this error is the same as the global error stored
208219
print('Average freq-by-freq error:\t {:1.3f}'.format(np.mean(errs_fm)))
209220
print('specparam model fit error: \t\t {:1.3f}'.format(fm.error_))

0 commit comments

Comments
 (0)