One option for analysing your long format data would be to use linear mixed effects modelling via the lme() function in the nlme package in R.
For example, a so-called random intercept linear mixed effects model for your data could be set up as follows:
library(nlme)
M1.lme <- lme(DV ~ IV1*IV2, random = ~ 1|S.ID, data = DATASET, method = "REML")
A linear mixed effects model with random intercepts and random slopes for IV2 can be specified as:
M2.lme <- lme(DV ~ IV1*IV2, random = ~ 1 + IV2|S.ID, data = DATASET, method = "REML")
You have to pay attention to what fitting method you will use for these models - REML (Restricted Maximum Likelihood) or ML (Maximum Likelihood). See, for example, How to decide whether to set REML to True or False?.
The fact that the variables DV1 and DV2 are allowed to interact in your model (e.g., Model M1) is specified via the term IV1*IV2. An equivalent but explicit way to specify the main effects of DV1 and DV2 and the interaction between DV1 and DV2 would be as follows:
M1.lme <- lme(DV ~ IV1 + IV2 + IV1:IV2, random = ~ 1|S.ID, data = DATASET, method = "REML")
For either model, you can construct an analysis of deviance table using these commands:
library(car)
Anova(M1)
Anova(M2)
You can also conduct post-hoc/contrast analyses, as explained for instance here: http://rcompanion.org/handbook/I_09.html.
The same types of models can also be fitted with the lmer() function from the lme4 package, though the syntax is slightly different (see https://www.r-bloggers.com/how-to-do-repeated-measures-anovas-in-r/):
library(lme4)
library(lmerTest)
M1.lmer <- lmer(DV ~ IV1*IV2 + (1|S.ID), data = DATASET, REML = TRUE)
anova(M1.lmer)
M2.lmer <- lmer(DV ~ IV1*IV2 + (1 + IV2|S.ID), data = DATASET, REML = TRUE)
anova(M1.lmer)
Note that the main effects of IV1 and IV2 and the interaction between them are assumed to be the same across all subjects in the model M1.lme (or M1.lmer). However, model M2.lme (or M2.lmer) assumes the main effect of DV2 to differ across subjects. In principle, it is also possible to allow the interaction between DV1 and DV2 to differ across subjects, though most people seem not to consider that possibility for some reason (likely related to the increased complexity of the model):
M3.lme <- lme(DV ~ IV1 + IV2 + IV1:IV2, random = ~ 1 + IV2 + IV1:IV2|S.ID, data = DATASET, REML = TRUE)
Anova(M1.lme)
M3.lmer <- lme(DV ~ IV1 + IV2 + IV1:IV2 + (1 + IV2 + IV1:IV2|S.ID), data = DATASET, REML = TRUE)
anova(M1.lme)
Anyway, this should give you enough to get started!