Identifying distinct occupations with pseudo-Bayes Factors

[Ben Marwick]

Hi everyone,

I'm trying to identify whether or not I have distinct occupations in a set of radiocarbon ages. I want to know if I should consider them a single phase or if there is a gap between them. I've come up with four possible models, and got A_model scores. The model scheme is meant to be like this table below, where 1-3 are phase names, and each row is a date (each column is a model), so model A is a single phase, model D has three phases, etc.:

A, B, C, D

1, 1, 1, 1

1, 1, 1, 1, 

1, 1, 2, 2, 

1, 2, 2, 3, 

1, 2, 2, 3, 

Could I ask for some feedback on if I have correctly implemented these models in the code below? And am I interpreting the A_model values correctly to indicate that model A is the best one in this case? (I took inspriation from this Q&A: https://groups.google.com/d/msg/oxcal/YRGxmMKWIIo/lOUbzBsP7_oJ). I would like to compute likelihood ratio between pairs of my models to compare with some published results. I've had a go at the pseudo-Bayes Factors using Andrew's instructions at https://groups.google.com/d/msg/oxcal/mDVhh8FNnGY/6ugk93lyFrEJ My attempt is below. I hope they're on the right track. I'd be most grateful for any guidance. 

Here are the models:

  ################## A ########### A=99

 

 Plot()

 {

  Sequence()

  {

   Boundary("Start 1");

   Phase("1")

   {

  C_Date("FC455", 187, 101);

  C_Date("FC486", 392, 58);

  C_Date("FC29", 822, 61);

  C_Date("FC20", 1789, 75);

  C_Date("FC460", 1916, 66);

   };

   Boundary("End 1");

  };

 };

 

 

  ################## B ########### A=97.1

 

 Plot()

 {

  Sequence()

  {

   Boundary("Start 1");

   Phase("1")

   {

    C_Date("FC455", 187, 101);

    C_Date("FC486", 392, 58);

    C_Date("FC29", 822, 61);

   };

   Boundary("trans");

   Phase("2")

   {

    C_Date("FC20", 1789, 75);

    C_Date("FC460", 1916, 66);

   };

   Boundary("End 1");

  };

 };

 

 ################## C ########### A=93.8

 Plot()

 {

  Sequence()

  {

   Boundary("Start 1");

   Phase("1")

   {

    C_Date("FC455", 187, 101);

    C_Date("FC486", 392, 58);

   };

   Boundary("trans");

   Phase("2")

   {

    C_Date("FC29", 822, 61);

    C_Date("FC20", 1789, 75);

    C_Date("FC460", 1916, 66);

   };

   Boundary("End 1");

  };

 };

 ################## D ########### A=90.2

  Plot()

 {

  Sequence()

  {

   Boundary("Start 1");

   Phase("1")

   {

    C_Date("FC455", 187, 101);

    C_Date("FC486", 392, 58);

   };

   Boundary("trans1");

   Phase("2")

   {

    C_Date("FC29", 822, 61);

   };

   Boundary("trans2");

   Phase("3")

   {

    C_Date("FC20", 1789, 75);

    C_Date("FC460", 1916, 66);

   };

   Boundary("End 1");

  };

 };

Here's where I compute the Bayes Factors

##### Pseudo-Bayes Factors  (using R here) ###

# From OxCal

A_A=99

B_A=97.1

C_A=93.8

D_A=90.2

# n= 6 dates for all models

F_model_A <- (A_A/100)^sqrt(6)

F_model_B <- (B_A/100)^sqrt(6)

F_model_C <- (C_A/100)^sqrt(6)

F_model_D <- (D_A/100)^sqrt(6)

# compute F_model as an approximation to a pseudo-Bayes factor

> (F_model_A / F_model_B) ^ sqrt(6)

[1] 1.1233

> (F_model_A / F_model_B) ^ sqrt(6)

[1] 1.1233

> (F_model_A / F_model_C) ^ sqrt(6)

[1] 1.382274

> (F_model_A / F_model_D) ^ sqrt(6)

[1] 1.748123

> (F_model_B / F_model_C) ^ sqrt(6)

[1] 1.230547

> (F_model_B / F_model_D) ^ sqrt(6)

[1] 1.556238

> (F_model_C / F_model_D) ^ sqrt(6)

[1] 1.264672

# Conclusion: Looks like A:D is the most extreme, gaps do exist in the dates.

Is that basically on the right track?

thanks very much,

Ben

-- 
Ben Marwick, Assistant Professor, Department of Anthropology
Box 353100, University of Washington
Seattle, WA 98195-3100 USA

(+1) 206.552.9450 bm...uw.edu @benmarwick
http://faculty.washington.edu/bmarwick/ 

[MILLARD A.R.]

Ben,

Your Pseudo-Bayes factors range only 1.1 to 1.7. As I noted in my original post, Bayes Factors greater than 5 or 10 are generally reckoned to be strong evidence to prefer one model over another. Jeffreys would place yours in the "barely worth mentioning" category.

Wikipedia has a good summary of the two major schemes for interpreting Bayes Factors. https://en.wikipedia.org/wiki/Bayes_factor#Interpretation


More generally, I'd still like to see some comment on the mathematical validity of my original proposal. It could also do with some testing on simulated datasets.




Best wishes

Andrew
--
 Dr. Andrew Millard 
e: A.R.M...@durham.ac.uk | t: +44 191 334 1147
 w: http://www.dur.ac.uk/archaeology/staff/?id=160
 Senior Lecturer in Archaeology, Durham University, UK

> --
> You received this message because you are subscribed to the Google
> Groups "OxCal" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to oxcal+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Rayfo...@aol.com]

Andrew,

 

Can you clarify for me why we are using A_model indices and '^sqrt(n)' in the Bayes factor calculations, or is it only to show the maths, where you say:

 

"Thus using that F_model = (A_model/100)^sqrt(n) where n is the number of dates, for two models A and B, we have that F_modelA/F_modelB =(A_modelA/AmodelB)^sqrt(n) and this can be used as an approximation to a Bayes Factor between the two models. Bayes Factors greater than 5 or 10 are generally reckoned to be strong evidence to prefer one model over another. "

 

What I'm getting at is that: 'View - Model Specification' gets to Fmodel directly viz:

 

 

Then in Excel

 

Bayes Factor Calcs			n- No of dates		40
			Sqrt (n)		6.324555
	F_modelA	3967		A_modelA	370.66
	F_modelB	52.6		A_modelB	187.12
	Bayes factor	75.41871		F_modelA/F_modelB=(A_modelA/A_modelB)^sqrt(n)

 

 

The above two models are a comparison of a sequence of 20 dates with (a) Lognormal intervals and (b) a P sequence (Poisson intervals).  So in this case, ceteris paribus, the model with lognormal intervals is supported.

 

I'm sorry, I just noticed your May requests for comments.......

 

Best wishes

 

Ray

 

 

In a message dated 10/07/2015 08:45:13 GMT Daylight Time, a.r.m...@durham.ac.uk writes:

https://en.wikipedia.org/wiki/Bayes_factor#Interpretation

[MILLARD A.R.]

Ray,

> Can you clarify for me why we are using A_model indices
> and '^sqrt(n)' in the Bayes factor calculations

...

> 'View - Model Specification' gets to Fmodel directly

I can’t remember why I gave that formula. I suspect the original question was about how to compare models using A_model, so that was my starting point. Or perhaps I forgot that F_model was directly available from OxCal.

[Rayfo...@aol.com]

Andrew,

 

Thanks for clearing that up.

 

I've been comparing models of a Sequence A with models of a Sequence B with inserted intervals (other than that allowed by the program).  However the 'n' doubles to accommodate the inserted Intervals so the Bayes Factor Fmodel comparison is fraught.  By inserting a nominal interval in Sequence A, I can get the same 'n' which allows me to compare alternative intervals between events.

 

Just toying.

 

Best wishes

 

Ray

 

[Ben Marwick]

Hi Andrew and Ray,

Thanks for taking a look and your comments, that's all very helpful.

best,

Ben

[MILLARD A.R.]

Yesterday a Google Scholar update led me to one comment on the idea of using F_model to compare models, which is worth noting. This is in a recent book chapter by Christopher:

https://books.google.co.uk/books?id=pi6sCQAAQBAJ&lpg=PA272&ots=WU0LjdRDqv&lr&pg=PA290#v=onepage&q=tempted&f=false

Best wishes

Andrew
--
 Dr. Andrew Millard 
e: A.R.M...@durham.ac.uk | t: +44 191 334 1147
 w: http://www.dur.ac.uk/archaeology/staff/?id=160
 Senior Lecturer in Archaeology, Durham University, UK

> -----Original Message-----
> From: ox...@googlegroups.com [mailto:ox...@googlegroups.com] On Behalf
> Of Ben Marwick

> Sent: 11 July 2015 08:35
> To: ox...@googlegroups.com
> Subject: Re: Identifying distinct occupations with pseudo-Bayes Factors
>
> Hi Andrew and Ray,
>
> Thanks for taking a look and your comments, that's all very helpful.
>
> best,
>
> Ben
>
> On Friday, July 10, 2015 at 1:06:46 PM UTC-7, Ray wrote:
>
>
> Andrew,
>
> Thanks for clearing that up.
>
> I've been comparing models of a Sequence A with models of a
> Sequence B with inserted intervals (other than that allowed by the
> program). However the 'n' doubles to accommodate the inserted Intervals
> so the Bayes Factor Fmodel comparison is fraught. By inserting a
> nominal interval in Sequence A, I can get the same 'n' which allows me
> to compare alternative intervals between events.
>
> Just toying.
>
> Best wishes
>
> Ray
>
> In a message dated 10/07/2015 15:05:38 GMT Daylight Time,

> a.r.m...@durham.ac.uk <javascript:> writes:
>
> Ray,
>
> > Can you clarify for me why we are using A_model indices
> > and '^sqrt(n)' in the Bayes factor calculations
> ...
> > 'View - Model Specification' gets to Fmodel directly
>
> I can’t remember why I gave that formula. I suspect the
> original question was about how to compare models using A_model, so that
> was my starting point. Or perhaps I forgot that F_model was directly
> available from OxCal.
>
>
>
> Best wishes
>
> Andrew
> --
> Dr. Andrew Millard

> e: A.R.M...@durham.ac.uk <javascript:> | t: +44 191 334

> 1147
> w: http://www.dur.ac.uk/archaeology/staff/?id=160
> Senior Lecturer in Archaeology, Durham University, UK
>
>
> --
> You received this message because you are subscribed to the
> Google Groups "OxCal" group.
> To unsubscribe from this group and stop receiving emails from

> it, send an email to oxcal+un...@googlegroups.com <javascript:> .

[Gianmarco Alberti]

Hello,

I was following this interesting thread, and I wanted to be sure to have got it right:

for the purpose of models comparison, one can use or the Fmodel value reported by OxCal (in this case the Bayesian Factor would be the ratio between the Fmodel values of two models; OR, one can use the Amodel reported in the OxCal output table (in this case the BF would be (Amodel_A/Amodel_B)/sqrt(nos of dates)), right?

Best 

Gm

Inviato da iPad

--

[Rayfo...@aol.com]

Ciao Gianmarco,

 

Stai bene?

 

Not quite.   'n' is the number of parameters in the model calculation (I think between the Start and End boundaries) this is not necessarily limited to the number of dates, there may be other parameters taken into account.

For instance this model has 8 dates with intervals so 'n' is 17.

 

 

Also the formula using Amodel index is the ratio of the two Amodel indices raised to the power of Sqrt(n). (not divided by Sqrt(n)).    Using Fmodel indices it is just the ratio of the two Fmodel indices.

 

However it is a pseudo Bayes factor and the point made by Andrew on the caveat in Christopher's paper is well taken.

 

Best wishes,

 

Ray

[Gianmarco Alberti]

Hello Ray,

Thanks for the clarification. The reference to the 'division' was an error on my part, sorry. I was meaning: ^sqrt(n).

Caveats noted!!!

Best

Gm

Inviato da iPad

[MILLARD A.R.]

The really important thing is that the null model part of the F_model is the same when comparing two models. That means that the same n parameters must all be in both models and have the same likelihoods. If there are parameters other than dates this will need some careful thought to make sure the two are comparable. For Ben’s models this is true.

 

Ray, what type of model are you constructing that has parameters other than dates?

 

Best wishes

Andrew
--
 Dr. Andrew Millard 

   e: A.R.M...@durham.ac.uk | t: +44 191 334 1147
   w: http://www.dur.ac.uk/archaeology/staff/?id=160

 Senior Lecturer in Archaeology, Durham University, UK

 

From: ox...@googlegroups.com [mailto:ox...@googlegroups.com]
Sent: 11 July 2015 20:32
To: ox...@googlegroups.com
Subject: Re: Identifying distinct occupations with pseudo-Bayes Factors

 

Ciao Gianmarco,

[Rayfo...@aol.com]

Hi Andrew,

 

I'm toying with a few special groupings but as an example of models with parameters other than only dates, then V_Sequence and V_Sequence Equivalent from the Help file is a fair example:

 

"The V_Sequence method, carried forward from previous versions of OxCal extends the D_Sequence methodology to allow for uncertainty in the gaps between events. The events are also constrained to be in order (so the Normal uncertainty is truncated at zero). In the current version the same can be achieved by applying a prior to the interval between events in a normal sequence. Both of the following are equivalent."

This results in the V_Sequence having n=8 and the V_Sequence Equivalent n=17.  The Gaps N(10,5)  do not appear as parameters, whereas in the V_Squence Equivalent the Intervals N(10,5),  do.

 

Below: Although the Amodel indices are the same, the Fmodel indices are not.    The Null model parameters are not the same so Pseudo Bayes factor won't be valid using the Amodel ratios, but it serves to illustrate that n is not only a count of the dates in the model.

 

However, since the 1/sqrt(n) appears in the calculation of Amodel from Fmodel, is the 'n' a necessary element in the Amodel ratio Bayes Factor calculation, but not in the Fmodel ratio calculation?  I certainly don't know.  I guess if I wanted to frame the question it would be something like:

 

Can I compare two Fmodel indices by their ratio, but if I want to use Amodel indices I need to know that n is the same in both models? 

 

 

V_Sequence V_Sequence() { Boundary("Start"); Gap(10,5); R_Date("A",2023,20); Gap(10,5); R_Date("B",1961,20); Gap(10,5); R_Date("C",1999,20); Gap(10,5); R_Date("D",1966,20); Gap(10,5); R_Date("E",1954,20); Gap(10,5); R_Date("F",1936,20); Gap(10,5); R_Date("G",1948,20); Gap(10,5); R_Date("H",1925,20); Gap(10,5); Boundary("End"); }; Best wishes   Ray

V_Sequence equivalent Sequence() { Boundary("Start"); Interval(N(10,5)); R_Date("A",2023,20); Interval(N(10,5)); R_Date("B",1961,20); Interval(N(10,5)); R_Date("C",1999,20); Interval(N(10,5)); R_Date("D",1966,20); Interval(N(10,5)); R_Date("E",1954,20); Interval(N(10,5)); R_Date("F",1936,20); Interval(N(10,5)); R_Date("G",1948,20); Interval(N(10,5)); R_Date("H",1925,20); Interval(N(10,5)); Boundary("End"); };

[MILLARD A.R.]

Hi Ray,

 

The important thing is that the likelihoods are the same in both models and only the priors change. So not only must there be the same number of parameters, but they must have the same likelihood functions so that the null model is the same. If you compare the list of likelihoods for the two v_sequence models you can see they are different so the  element of the F_model is not the same and does not cancel when you take the ratio of the two F_model values. If I add additional boundaries to the second model, the likelihood function does not change (though the labelling of the parameters does) and the F_model values can be compared.

 

(very wide table follows)

 

V_Sequence	V_Sequence Equivalent	V_Sequence Equivalent modified

V_Sequence()   {    Boundary("Start");    Gap(10,5);    R_Date("A",2023,20);    Gap(10,5);    R_Date("B",1961,20);    Gap(10,5);    R_Date("C",1999,20);    Gap(10,5);    R_Date("D",1966,20);    Gap(10,5);    R_Date("E",1954,20);    Gap(10,5);    R_Date("F",1936,20);    Gap(10,5);    R_Date("G",1948,20);    Gap(10,5);    R_Date("H",1925,20);    Gap(10,5);    Boundary("End");   };

Sequence()

Boundary("Start pause");    Boundary("End pause");

Interval(N(10,5));    R_Date("E",1954,20);    Interval(N(10,5));    R_Date("F",1936,20);    Interval(N(10,5));    R_Date("G",1948,20);    Interval(N(10,5));    R_Date("H",1925,20);    Interval(N(10,5));    Boundary("End");   };

[Rayfo...@aol.com]

Hi Andrew,

 

Thanks for taking the time over this.  I think I follow.  However I am puzzled in that the V_Seq and V_Seq_Equiv have the same Amodel (about 104) when Gap N(10,5) and Interval N(10,5) are used.  But this does not hold when I use another Gap or Interval, say Gap N(12,6) and Interval N(12,6).  So there is something else going on that I don't follow.  Can the Amodel agreement be somewhat fortuitous in the Help example?

 

Best wishes

 

Ray

[MILLARD A.R.]

> From: ox...@googlegroups.com [mailto:ox...@googlegroups.com]
> Sent: 15 July 2015 20:35

>
> Thanks for taking the time over this. I think I follow. However I am
> puzzled in that the V_Seq and V_Seq_Equiv have the same Amodel (about
> 104) when Gap N(10,5) and Interval N(10,5) are used. But this does not
> hold when I use another Gap or Interval, say Gap N(12,6) and Interval
> N(12,6). So there is something else going on that I don't follow. Can
> the Amodel agreement be somewhat fortuitous in the Help example?

Ray,

I think it must just be a chance match. When I ran them they weren't exactly the same.

[Rayfo...@aol.com]

Andrew,

 

Thanks for all this.

 

Can I just summarise so that my understanding is correct?

 

In OxCal, two Fmodel indices may be compared provided:

 

1. 'n'  the number of likelihoods is the same in both models.

 

2. Each individual Likelihood (the pdf) is identical in both models.

 

3.  These identical Likelihoods are, in addition to the R_Dates, C_Dates and Dates, also all other Likelihoods such as created by Simulated dates, Intervals with expression, Span with expression, in fact any parameter that has a valid expression included that generates a Likelihood pdf.

 

4.  Only the priors in the models may be changed if the Fmodel comparisons are to remain valid.

 

5. The magnitude of the Fmodel ratio has a bearing on its significance (small s)

 

Is it worth emphasizing that any two models (even with different priors) that include Simulated dates should not be compared using Fmodel ratios since each run of a model with R_Simulate, etc will generate a different Likelihood?

 

On a different, but similar vein, can you comment on the following:

 

Assuming a model with a priori interval pdfs between the dates.  The model is run umpteen times with only the Interval Mean changed.  Then the Fmodel and Amodel vs the Interval Mean is plotted.  The Likelihoods change with each new Interval, as do the Interval posteriors, so the Fmodel index can't be compared, however the plot shows a distinct Fmodel and Amodel increase around certain Interval values.

 

Can anything be said regarding which range of Interval is preferred? (say 14 to 18, the posterior will be less)

 

Interval (LnN(ln(x),ln(2)))
x=	Fmodel	Amodel
6	1.8498	116.09
8	2.47702	124.61
10	4.53871	144.32
12	6.59131	157.99
14	8.22518	166.71
16	9.05746	170.65
18	9.0458	170.6
20	8.10027	166.09
22	7.0429	160.55
24	5.46262	150.95

 

 

Best wishes

 

Ray