Reading the files

To run the example, download the 3 files in https://github.com/danilinares/quickpsy/tree/master/inst/extdata/example1 and add them to your working directory.

These files have a quite standard tidy form in which each column corresponds to a variable and each row is an observation.

The files can be read manually using the R base read.table function

library(quickpsy)

dat1 <- read.table('aa1.txt', head = T)
dat2 <- read.table('bb1.txt', head = T)
dat3 <- read.table('cc2.txt', head = T)

dat1 <- dat1 %>% mutate(subject = 'aa')
dat2 <- dat2 %>% mutate(subject = 'bb')
dat3 <- dat3 %>% mutate(subject = 'cc')

dat <- bind_rows(dat1, dat2, dat3)

but quickpsy facilitates its reading with the function quickreadfiles

library(quickpsy)
dat <- quickreadfiles(subject = c('aa', 'bb', 'cc'),
                      session = c('1', '2'))

Fitting

We use the main function of quickpsy, which has the name of the package, to do the fit. We introduce the data, the name of the explanatory variable, the name of the response variable and the name of the grouping variables.

fit <- quickpsy(dat, phase, resp, grouping=.(ecc, subject)) 
# we can also introduce the name of the arguments: 
# fit <- quickpsy(d = dat, x = phase, k = resp, grouping=.(ecc, subject))

quickpsy creates a list of class quickpsy with all the results. We can use names to obtain the names of the output variables

names(fit)
##  [1] "x"                    "k"                    "n"                   
##  [4] "guess"                "lapses"               "averages"            
##  [7] "groups"               "funname"              "log"                 
## [10] "psyfunguesslapses"    "limits"               "pariniset"           
## [13] "parini"               "optimization"         "par"                 
## [16] "pariniset"            "ypred"                "curves"              
## [19] "sse"                  "thresholds"           "logliks"             
## [22] "loglikssaturated"     "deviance"             "avbootstrap"         
## [25] "parbootstrap"         "logliksboot"          "logliksbootsaturated"
## [28] "devianceboot"         "aic"                  "parcomparisons"      
## [31] "curvesbootstrap"      "thresholdsbootstrap"  "thresholdcomparisons"

To obtain the parameters, for example, we can do

parameters <- fit$par
parameters
## Source: local data frame [18 x 6]
## Groups: ecc, subject [?]
## 
##      ecc subject   parn       par    parinf    parsup
##    <int>  <fctr> <fctr>     <dbl>     <dbl>     <dbl>
## 1      1      aa     p1 195.70563 181.82578 207.84797
## 2      1      aa     p2 124.18318 110.15372 139.23549
## 3      1      bb     p1 216.78897 206.29337 227.53820
## 4      1      bb     p2  71.65806  61.78054  81.55831
## 5      1      cc     p1 136.09043 121.49481 152.10270
## 6      1      cc     p2 178.92481 161.63406 208.42624
## 7      2      aa     p1 122.68393 112.21748 133.43589
## 8      2      aa     p2  79.93172  72.25307  89.03789
## 9      2      bb     p1 149.89178 141.63102 160.02312
## 10     2      bb     p2  54.49609  46.90926  63.06072
## 11     2      cc     p1  53.77198  38.77197  66.84737
## 12     2      cc     p2 114.46046  99.50036 132.54332
## 13     3      aa     p1  72.07884  62.98318  82.26360
## 14     3      aa     p2  67.12062  57.33812  75.59393
## 15     3      bb     p1  95.34767  83.93705 104.58408
## 16     3      bb     p2  56.54547  46.92185  66.15872
## 17     3      cc     p1  77.03154  65.65874  90.59112
## 18     3      cc     p2  98.64317  85.31828 110.82372

Plotting

To plot the curves, for example, we use the quickpsy function plotcurves

plot1 <- plot(fit)
plot1

As plot1 is a ggplot, we can modify it accordingly. For example

plot1 + theme_classic() + ylab('Proportion') + theme(panel.border = element_blank())

We can also use the data frames in the output of quickpsy to make the plot using ggplot2 commands directly

 ggplot() +
  facet_wrap(~subject) +
  geom_line(data = fit$curvesbootstrap, aes(x = x, y = y, color = factor(ecc), 
                    group=paste(sample,ecc)), lwd = .2, alpha = .4) +   
  geom_line(data = fit$curves, aes(x = x, y = y, group = factor(ecc)), 
            color = 'black') +
  geom_point(data = fit$averages, aes(x = phase, y = prob, color = factor(ecc)), size = 4)

Finally, we plot the parameters and the tresholds (plotting different ecc in different colors for the thresholds)

plotpar(fit)
## Warning: Ignoring unknown aesthetics: fill

plotthresholds(fit, color = ecc)
## Warning: Ignoring unknown aesthetics: fill