Latent Learning of Skill in Dance Dance Revolution I'm trying to model and plot the rate of learning Dance Dance Revolution over the past year.
I have data of the format below, with the score received for a song on a specific date and difficulty level. Songs have been played many times over the year, at a variety of difficulty levels. I'd like to model the skill learned over the year.




Date Played
Song
Difficulty Level
Game Version
Score
Letter Grade
Num Perfects
Num Greats
Num Goods
Num Almosts
Num Boos
Num O.K.s
Max Combo




March 12, 2021
Girls Just Wanna Have Fun
Expert
DDR Supernova
8381742
A
153
52
4
0
7
23
69




Each song is unique per game version, and the scores across game versions are non standard (some are in the hundred thousand range and some in the thousand range). There are a total of 8 game versions, 5 difficulty levels, and ~100 songs.
A trivial solution would be to do a linear regression of the score over days, but this would not take into account differences by song, difficulty level, or game version.
Any ideas how to model the learning rate in a less trivial way?
 A: I suggest a simple generalized linear model to start.
$score = \beta_0 + \beta_1 days + \beta_2 song + \beta_3 difficulty + \beta_4 version$
where

*

*day = days since start of training

*song = a categorical encoding of the song

*difficulty = a categorical encoding of the difficulty

*version = a categorical encoding of the version

Since score is positive, and since the residuals might be long tailed to higher scores, I chose a lognormal error distribution.  You will need to inspect the data to determine the right error distribution.
I didn't add player to the regression, but you could add a categorical for player also.
in R:

require(ggplot2)

dates <- as.Date("2021-01-01") + 0:365
song <- c(LETTERS, letters)
difficulty <- paste0("level", 1:5)
version <- paste0("version", 1:8)

# simulate data, you would read your data in (read.csv or something like it)
set.seed(1934493)
N <- 300
dat <- data.frame(date = dates[sample(1:length(dates), size = N, replace = TRUE)],
                  song = factor(song[sample(1:length(song), size = N, replace = TRUE)]),
                  diff = factor(difficulty[sample(1:length(difficulty), size = N, replace = TRUE)]),
                  vers = factor(version[sample(1:length(version), size = N, replace = TRUE)]))
dat$dayssincestart <- as.numeric(dat$date - as.Date("2021-01-01"))
dat$score <- with(dat, dayssincestart * 10 + as.numeric(song) * 1 + as.numeric(diff) * 2 + as.numeric(vers) * 3 + rlnorm(N, 1, 1))
# ensure that each song, diff, and vers is done at least 2x
all(table(dat$song) >= 2)
all(table(dat$diff) >= 2)
all(table(dat$vers) >= 2)

glm1 <- glm(score ~ dayssincestart + song + diff + vers, data = dat, family = gaussian(link = "log"))
glm0 <- glm(score ~ 1, data = dat, family = gaussian(link = "log"))

# accounting for song, difficulty, and version, is there a score trend with time?
#   (1) look at the model to see it it is significant overall
anova(glm1, glm0, test = "LRT")
# Yes, P < 0.05, there is at least one significant relationship
#   (2) look at the coefficient on days_since_start to determine if you learning with time
summary(glm1)
# Yes, p < 0.05 and the estimate is positive (5.29E-3 increase in score per day after accounting for level, song, and version)

# plot data
ggplot(dat, aes(x = date, y = score, col = diff)) +
  geom_point() +
  labs(x = "Date", y = "Score", col = "Difficulty")

# plot diagnostics
plot(glm1, which = 1) # structure of the residuals indicates a problem which would require re-fitting
plot(glm1, which = 2)
plot(glm1, which = 3)
plot(glm1, which = 4)


