The same result (given in Coolwater's answer) can be achieved using NMinimize without initial parameter guess:

sol = ParametricNDSolveValue[{ S'[t] == -a*j[t]*S[t]/12500, q'[t] == a*j[t]*S[t]/12500 - 0.2*q[t], j'[t] == 0.2*q[t] - b*j[t], r'[t] == b*j[t], F'[t] == 0.2*q[t], S[0] == 12500, q[0] == 1, j[0] == 10, r[0] == 0, F[0] == 0}, F, {t, 0, 338}, {a, b}, DependentVariables -> {S, q, j, r, F}]; 

NMinimize shows several warnings

mini = NMinimize[Total@Map[(sol[a, b][#[[1]] ] - #[[2]])^2 &, data],{a, b}] (* {766228., {a -> 0.220993, b -> 0.186785}} *) 

but the optimal parameter pair is the same!

Show[{Plot[ sol[a, b][t] /. mini[[2]] , {t, 0, 338}], ListPlot[data]}] 

remark:

FindFit only works with Method->"NMinimize"

fit = Quiet@FindFit[data, sol[a, b][t], {a, b}, t, Method -> "NMinimize"] 

The same result (as given in Coolwater's answer) can be achieved using NMinimize without initial parameter guess:

sol = ParametricNDSolveValue[{ S'[t] == -a*j[t]*S[t]/12500, q'[t] == a*j[t]*S[t]/12500 - 0.2*q[t], j'[t] == 0.2*q[t] - b*j[t], r'[t] == b*j[t], F'[t] == 0.2*q[t], S[0] == 12500, q[0] == 1, j[0] == 10, r[0] == 0, F[0] == 0}, F, {t, 0, 338}, {a, b}, DependentVariables -> {S, q, j, r, F}]; 

NMinimize shows several warnings

mini = NMinimize[Total@Map[(sol[a, b][#[[1]] ] - #[[2]])^2 &, data],{a, b}] (* {766228., {a -> 0.220993, b -> 0.186785}} *) 

but the optimal parameter pair is the same!

Show[{Plot[ sol[a, b][t] /. mini[[2]] , {t, 0, 338}], ListPlot[data]}] 

remark:

FindFit only works with Method->"NMinimize"

fit = Quiet@FindFit[data, sol[a, b][t], {a, b}, t, Method -> "NMinimize"] 

The same result (given in Coolwater's answer) can be achieved using NMinimize without initial parameter guess:

sol = ParametricNDSolveValue[{ S'[t] == -a*j[t]*S[t]/12500, q'[t] == a*j[t]*S[t]/12500 - 0.2*q[t], j'[t] == 0.2*q[t] - b*j[t], r'[t] == b*j[t], F'[t] == 0.2*q[t], S[0] == 12500, q[0] == 1, j[0] == 10, r[0] == 0, F[0] == 0}, F, {t, 0, 338}, {a, b}, DependentVariables -> {S, q, j, r, F}]; 

NMinimize shows several warnings

mini = NMinimize[Total@Map[(sol[a, b][#[[1]] ] - #[[2]])^2 &, data],{a, b}] (* {766228., {a -> 0.220993, b -> 0.186785}} *) 

but the optimal parameter pair is the same!

Show[{Plot[ sol[a, b][t] /. mini[[2]] , {t, 0, 338}], ListPlot[data]}] 

The same result (given in Coolwater's answer) can be achieved using NMinimize without initial parameter guess:

sol = ParametricNDSolveValue[{ S'[t] == -a*j[t]*S[t]/12500, q'[t] == a*j[t]*S[t]/12500 - 0.2*q[t], j'[t] == 0.2*q[t] - b*j[t], r'[t] == b*j[t], F'[t] == 0.2*q[t], S[0] == 12500, q[0] == 1, j[0] == 10, r[0] == 0, F[0] == 0}, F, {t, 0, 338}, {a, b}, DependentVariables -> {S, q, j, r, F}]; 

NMinimize shows several warnings

mini = NMinimize[Total@Map[(sol[a, b][#[[1]] ] - #[[2]])^2 &, data],{a, b}] (* {766228., {a -> 0.220993, b -> 0.186785}} *) 

but the optimal parameter pair is the same!

Show[{Plot[ sol[a, b][t] /. mini[[2]] , {t, 0, 338}], ListPlot[data]}] 

remark:

FindFit only works with Method->"NMinimize"

fit = Quiet@FindFit[data, sol[a, b][t], {a, b}, t, Method -> "NMinimize"]