I usually use a pseudo-object-oriented paradigm in which an object has the form classHead[data]. Methods are defined for classHead as subvalues. I usually use an Association[] for data, but pre V10 this was not possible. Here are two variants. If I want to update data in my object (the mutable variant), then I need to be able to modify the association. To do this the argument has to be a symbol that refers to an association. In a package, I often create this symbol in a private subcontext, but below I do it manually in the "Global`" context. Further the argument must remain a symbol, so classHead[] must have an attribute like HoldAll. In the immutable variant, the data does not change, and the argument need not be held. In both variants, the methods can either return parts of the data or compute something from the data.

myObjImmutable // ClearAll;
myObjImmutable[data_Association]["Total"] := (* the "Total" method *)
  Total@myObjImmutable[data]["Data"];
myObjImmutable[data_Association][method_String] := (* generic data lookup *)
  Lookup[data, method, Message[myObjImmutable::notmethod]];

myObjMutable // ClearAll;
myObjMutable // Attributes = {HoldAll};
myObjMutable[data_Symbol]["Total"] /; AssociationQ@data :=  (* the "Total" method *)
  Lookup[data, "Total", data["Total"] = Total[data["Data"]]];
myObjMutable[data_Symbol][method_String] /; AssociationQ@data := (* generic data lookup *)
  Lookup[data, method, Message[myObjMutable::notmethod]];

my1 = myObjImmutable[<|"Data" -> Range[5]|>]
(*  myObjImmutable[<|"Data" -> {1, 2, 3, 4, 5}|>]  *)

my1["Total"]
(*  15  *)

my1 (* unchanged *)
(*  myObjImmutable[<|"Data" -> {1, 2, 3, 4, 5}|>]  *)

myData = <|"Data" -> Range[5]|>;
my2 = myObjMutable[myData]
(*  myObjMutable[myData]  *)

my2["Total"]
(*  15  *)

myData (* has been updated with "Total" *)
(*  <|"Data" -> {1, 2, 3, 4, 5}, "Total" -> 15|>  *)

So the method employed by LinearModelFit is a hybrid: It does return a function, FittedModel, but this is an existing function defined externally (and hidden) and is returned with a bunch of info bundled into an Association as its argument:

FittedModel[<|"Type" -> "Linear", 
 "Model" -> <|"FittedParameters" -> {0.186441, 0.694915}, 
   "IndependentVariables" -> {x}, "BasisFunctions" -> {1, x}, 
   "LinearOffset" -> <|"Function" -> 0, "Values" -> 0|>|>, 
 "KeyNames" -> Automatic, "Weights" -> <|"ExampleWeights" -> 1.|>, 
 "InputData" -> {{0, 1}, {1, 0}, {3, 2}, {5, 4}}, 
 "UserDefinedDesignMatrixQ" -> False, 
 "DesignMatrix" -> {{1., 0.}, {1., 1.}, {1., 3.}, {1., 5.}}, 
 "Localizer" -> 
  Function[Null, FittedModels`LocalizedBlock[{x}, #1], {HoldAll}], 
 "Options" -> {}|>]

This bundled info is enough for FittedModel to calculate any of the requested outputs. So FittedModel must be defined using SubValues something like this:

FittedModel[params_]["PValue"]:=...
FittedModel[params_]["ParameterTable"]:=

So I think I am beginning to see how it works! Don't understand the "Localizer" association yet, but little steps...

Ah — I look up SubValues and see what you mean:

myDataObj[x_, y_]["Properties"] := {"Properties", "Sum", "Product"}
myDataObj[x_, y_]["Sum"] := x + y
myDataObj[x_, y_]["Product"] := x y

and in fact the use of myDataObj as a name strongly suggest this is the 'official' way! But these functions are defined as having inputs: you can't just call myDataObj["Sum"], you have call myDataObj[x,y]["Sum"]. So in the typical use case the necessary pieces need to be defined like this

myDataObj["Product"] := x y

Here is how I got it to work:

fittingFunction[x_, y_] := (fitObject["Sum"] := x + y; fitObject["Product"] := x*y; fitObject)
myFit = fittingFunction[2, 3];
myFit["Sum"] (* Yields 5 *)
myFit["Product"] (* Yields 6 *)

Can anyone from Wolfram confirm that this is the general framework?

I note that one can usually also do Normal[myFit]; still need to figure that out!

Yes, I discovered this myself; should have done so earlier. The difference compared to my suggested approach is the FittedModel function is independently defined, not within the fitting function. The LinearModelFit function essentially does a 'first step' of calculations, and then passes a set of results from that (not necessarily corresponding to the possible queries, as you point out) to the second function as an association. The second function is presumably mostly a big Which[] based on all the possible string arguments.

Having figured this out, v15 is introducing a new 'unified' model fitting framework, so before implementing anything myself I guess I should grab the beta and see if anything significant has changed in how this part is handled!

An association works about this way:

myFit = Association["ParameterTable" -> whatever, 
   "AIC" -> somethingElse];
myFit["ParameterTable"]
myFit["AIC"]