Special Goods Production in UWO

Summary

Here I demonstrate that special goods production in UWO stands for: great success that produces a better new item than normal success does.

This can be a differently named item (such as master’s cannons) or the same item as normal success, but with better stats. It does not apply to improving (forging) existing items, to alchemy, or to making consumables or trade goods.

That this is so can be verified by crafting at sea with the Secrets of Craft sailor equipment (SE); it displays a message that its durability has decreased whenever such a special production is obtained, and never with other kinds of great success or any normal success. The durability loss on the SE increases with the required skill rank of the recipe, starting at 1 with 1 additional loss per 2 ranks (10 dura loss at r20).

I performed a factorial experiment to test the effects of special production EX and SE, as well as of casting rank and refining casting, by making 5,200 Minion (4) cannons, where the master’s version obtained with great success is “special goods”.

The results show that the effect of the SE is indeed to increase the degree (not rate) of great success. However, only the durability was increased, not the more important piercing power of the cannons. This makes the effect rather trivial and of limited usefulness. The effect can also be very small, but this depends on whether the casting skill is refined and on the casting rank. With refined r20 casting the effect was only around +2% durability, but was much larger at low casting rank and/or with unrefined casting (up to around +16%). The magnitude of the effect was similar to that of refining the skill and larger than increasing unrefined casting rank from r3 to r20. However, higher rank recipes may show different effect size than this r2 recipe did.

For piercing power, not only the SE but also refining the casting skill or increasing the casting rank from r3 to r20 had no effect. This was unexpected.

The effect of the special production EX on the other hand is to increase the rate (not degree) of great success. This effect was small, estimated at an additional 3% chance of a great success with Golden Arc EX r6 on.

Given the limited use cases for these items for special goods production, the high costs (dura loss), and the meagre effects, they are rather disappointing.

This is part 4 in a series of investigations about crafting in UWO.

Part 1 concerned alchemy experiments, part 2 making weapons with casting, and part 3 recipes for forging:

Stinker (2023) Report from a FABtory

Stinker (2024) Crafting success in UWO — experiments with Gorz

Stinker & Marsali (2024) Forging success in UWO — experiments with fur boots

What is special goods production?

In Uncharted Waters Online (UWO), the special goods production EX (Fig. 1) has no measurable effect in alchemy experiments (Stinker 2023), ordinary crafting of equipment (Stinker 2024), or forging of equipment (Stinker & Marsali 2024). This prompted a search for what “special goods” actually refers to.

FIGURE 1: The Golden Arc, special production EX r6

Instead of conducting further futile attempts to test when the EX does anything (and what it does), we have suggested (Stinker & Marsali 2024) that such an investigation should best be done by observing when the special production Sailor Equipment (SE, Fig. 2) loses durability, as this should only happen when the effect of the SE triggers (apart from the normal decay at sea over time). Further, we suggested that special goods production perhaps refers to when great success results in a qualitatively different item with a different name. Such items are produced from a limited number of recipes for ship parts, such as cannons, studding sails, corvus, etc.

FIGURE 2: The Secrets of Craft, r5 special production sailor equipment

I therefore explored many different recipes at sea with the sailor equipment active, to see when it loses durability (Tab. 1).

TABLE 1: Examples of production that did or did not result in durability loss of the Secrets of Craft sailor equipment when a great success was obtained, revealing what special goods production is and is not.

Type	Recipe		Product		Sailor Equip.
	Skill	Rank	Normal success	Great success1	Message	Loss2
special goods
equipment	casting	6	main gauche	(better stats)	yes	3
equipment	sewing	7	velvet ribbon	(better stats)	yes	3
equipment	casting	11	close helmet	(better stats)	yes	5
shippart	handicrafts	1	cedar plating	light cedar plating	yes	1
shippart	casting	2	minion (4)	master’s minion	yes	1
shippart	casting	3	corvus	(better stats)	yes	1
shippart	sewing	3	mizzen staysail	main staysail	yes	1
shippart	handicrafts	4	hardened red pine	(better stats)	yes	2
shippart	casting	5	large corvus	(better stats)	yes	2
shippart	casting	8	iron plating	light iron plating	yes	4
shippart	casting	10	light iron plating	wheatered iron plt.	yes	5
shippart	casting	13	hardened iron plt.	rolled iron plating	yes	6
shippart	casting	18	caprice (18)	master’s caprice	yes	9
not-so-special-goods
consumable	casting	2	spare rudder	(quantity)	no	—
consumable	handicrafts	3	lifesaver	(quantity)	no	—
forging	handicrafts	11	fur boots	(better stats)	no	—
trade goods	casting	5	cannon shot	(quantity)	no	—
trade goods	sewing	6	embroidery thread	(quantity)	no	—
trade goods	handicrafts	14	crystal	diamond	no	—
1 Denotes if great success resulted in a differently named product than normal success, or same product with better stats, or in larger quantity
2 Durability loss. This is per great success, and comes in addition to the normal loss over time at sea

Sailor equipment durability is indeed reduced at sea when a great success results in a qualitatively different item (Tab. 1).¹ There is even a message stating this when it happens (Fig. 3).

¹ Using sailor equip. maint. books against normal decay will not avoid this durability loss.

This message does not appear with normal success, nor with great success that produce consumables or trade goods or when forging an existing item, but invariably with great success when producing items (ship equipment or personal equipment). The great success does not need to result in a differently named item than from normal success; many recipes will instead produce the same item but with better stats. These will count as “special production” as well. Obtaining different trade goods from great success than from normal success (such as diamonds instead of crystals from white ore) does not count.

FIGURE 3: The message.

The message appears both when an item is produced from trade goods alone (such as a close helmet), or from upgrading an existing item into a better one (such as a long sword into a broadsword, or an iron plate into a rolled iron plate), but not when improving (i.e. forging) an existing item (e.g. fur boots).

Hence, it seems clear that special goods production means: great success that produces a new item with better stats than a normal success.

The amount of dura loss is not fixed, but is larger for higher rank recipes (Fig. 4).

FIGURE 4: The message from a rank 18 recipe.

Fig. 5 shows how the amount of dura loss per production is related to the required skill rank for the recipe. The observed values fits with an increased loss of 1 durability per 2 skill ranks.

FIGURE 5: The loss of durability (per production) of the special goods production Sailor Equipment increases with the required skill rank of the recipe. The observations (blue points) fits with an increased loss of 1 durability per 2 skill ranks, as indicated by the line in the figure.

What is the effect of the Special Goods Production EX and SE?

Now that we have clarified what special goods production is, we can ask what the EX and SE actually do.

Based on the in-game descriptions, the expected effect of the sailor equipment (Fig. 2) should be to increase the parameter value(s) of such “special products” (what we may call degree of success, or great success tier), and the effect of the EX (Fig. 1) to increase the frequency of getting “special goods” (i.e. the frequency of great success for such recipes, or what we may call rate of success).

Since getting as high stats as possible is often a goal when producing items, it would be of interest to know how large the effect of the sailor equipment is, for example compared to having max rank or refining the production skill.

It would also of course be of interest to increase the rate of great success, if that is what the EX does. However, an earlier study did not find a measurable effect of the EX on success rate when making the weapon Gorz, so the effect is likely small. Also, one can use a Master’s Secrets item to guarantee great success when doing expensive recipes.

Great success cannons (usually re-named “Master’s”) come in different tiers of great success for piercing power (PP) and the EX and Sailor Equipment might help getting high tier cannons.² Cannons should also be ideal for testing the effects of the SE and EX since they can be easy to produce in large quantity.

² Although the best cannons are from fixed recipes with NPCs at land, where the sailor equipment would have no effect (according to the description) and the EX would be irrelevant because one would use master’s secrets. The very best don’t even have a great success possibility and therefore no special goods production (or use of master’s secret) is possible.

I therefore designed a factorial experiment and made cannons.

For the rate of success (i.e. the EX) I was mostly concerned with detecting an effect at all (given the earlier failures to do so). For the degree of success (i.e. the SE) I was also interested in how large the effect is relative to the effects of skill rank and of refining the crafting skill.

Methods

The minion(4) recipe (Fig. 6) was chosen because it is the easiest cannon to make. Also, it minimizes the durability loss of the SE because of the low required rank (r2 casting), and maximizes the potential effect of having the maximum crafting rank vs. the required rank.

FIGURE 6: The recipe for crafting Minion (4) cannons

Production was made with a standardized set-up similar to previous reports in this series cited above. That is: always 4 Oxford production success skills on, Oxford item improvement on, aide in paymaster job and 120SS traits, no tarot active, Meister smith title on and apprentice nearby on land, and at sea in ship with workshop and the required management skill rank. Most production were done in Blacksmith job with Boston skill on to save mats, except when this gave too high casting rank.

In a factorial experiment I manipulated combinations of:

• refined or unrefined casting skill
• low (r3) or max (r20) casting rank
• with out without sailor equipment r5
• with or without EX r6

The experimental design was not complete for EX, in that no data was gathered for refined r3 casting with EX, or for unrefined r20 casting without EX. The remaining 12 combinations were all tested (Tab. 2).

For each cannon I recorded if the outcome was a normal or a great success (Fig. 7) and, in the case of great success, the piercing power and the durability.

FIGURE 7: Normal success and max great success minion (4).

TABLE 2: Treatment groups

Treatment	Where	Factors				Outcome		Extra2	Total
		Sailor Equip	Refined	Casting	EX	Normal	GS1
A	land	—	yes	r20	r6	142	118	1040	1300
B	land	—	yes	r20	no	138	122	1040	1300
C	land	—	no	r3	r6	115	145	0	260
D	land	—	no	r3	no	138	122	0	260
E	sea	r5	yes	r20	r6	125	135	0	260
F	sea	r5	yes	r20	no	141	119	0	260
G	sea	r5	no	r3	r6	129	131	0	260
H	land	—	no	r20	r6	140	120	0	260
I	land	—	yes	r3	no	130	130	0	260
J	sea	r5	yes	r3	no	134	126	0	260
K	sea	r5	no	r20	r6	113	147	0	260
L	sea	r5	no	r3	no	127	133	0	260
Total	—	—	—	—	—	1572	1548	2080	5200
1 GS = Great Success. This also gives the sample sizes for piercing power and durability.
2 Additional sample only used for success rate with and without EX

I made 260 cannons in each of the 12 treatment groups, in 13 batches of 20 cannons each using single production without closing the production window. Since the number of great success could vary between treatment groups, this means also the sample size for degree of success (PP and durability) vary somewhat.

FIGURE 8: Lowering the casting rank to r3 by mentoring

³ In mineral trader job to avoid extra ranks from blacksmith but retaining Meister title, and mentoring a trade level 0 character.

Even if the recipe is only r2, r3 was chosen for “low casting rank” since the lowest I could get with refined casting in favoured job³ was r15(-12), using mentoring to lower the rank (Fig. 8). To preserve identical conditions, I made sure the skill rank was still r3 with unrefined skill.

In addition I made some extra batches on land where only the outcome (normal or great success) was recorded. This was done to stabilize the estimates of success rate, lower the uncertainty, and increase the sample size and thus the statistical power to detect small differences in success rate. For these additional batches also 20 cannons were made in each batch, but with the “set amount” option. I did 52 such batches with and 52 without EX (on land, refined r20 casting).

In total, I obtained data on 5200 minion (4) cannons for success rate, and 1548 great success (out of 3120) gave data on piercing power and durability (Tab. 2).

I did not explicitly test the effect of workshop separate from sailor equipment (i.e. they were completely confounded since sailor equip requires workshop and I did not make any production at sea without sailor equipment). I make the assumption that any effect of producing at sea on success rate is from workshop and not from sailor equipment. Conversely, for degree of success I assume the effect is from sailor equipment and not from workshop.⁴

⁴ Three things change when producing at sea compared to land: 1) can no longer get the effect from Meister title with apprentice nearby on land (misses out on success bonus), but 2) can gain (a larger) success bonus from workshop, and 3) can get the effect of sailor equipment (needs workshop).

Statistical methods

I present the results in two ways: First visualizations of the raw data, with means and confidence intervals around the means and, for degree of success, also histograms. For rate of success, 95% binomial confidence intervals were obtained using the Wilson method. For degree of success, 95% confidence intervals were obtained by bootstrap resampling⁵ (adjusted (bias-corrected and accelerated) bootstrap percentile (BCa) method, 10K replicates).

⁵ https://en.wikipedia.org/wiki/Bootstrapping_(statistics)

⁶ https://en.wikipedia.org/wiki/Akaike_information_criterion

Second, I present a summary of statistical modelling of the data. Several candidate models were compared, and one selected (based on AIC⁶), and I visualize the predictions from the selected model.

Variation in success rate was modelled using logistic regression (see Note 1 below for an introduction).

Variation in piercing power and durability was modelled using linear models (function lm() in R). Since the residual variance was very heterogeneous I used robust (HC3) methods for the variance-covariance matrix to estimate uncertainty of the regression coefficients.⁷

⁷ The error structure could alternatively have been explicitly modelled, for example using generalized least squares, function gls() in the nlme library in R.

In Note 2 in the Appendix can be found some analyses of variation within batches to test if the RNG algorithm behaves as it should.

This report was written using Quarto. R code used to generate the results can be found in the Appendix together with software versions used. The raw data can be downloaded from the Appendix. Dark mode can be turned on/off in the top right, and figures can be enlarged by clicking on them.

Note 1: About odds, odds ratios and logistic regression

Most people (except gamblers) find it easier to think in terms of probability instead of odds, but the two concepts are just different scales of expressing the same thing.

Probability, \(p\), is an event’s fraction of all possible outcomes. That is, the expected proportion of times that an event will happen if repeated many times. Say, picking at random a card of hearts from a full standard 52-card deck will have \(p = 13/52 = 1/4 = 0.25\) (you can also express it as 1 in 4, or 25% if you really insist). Then \(1-p\) (often called \(q\)) will give the proportion of times that this event will not happen (i.e. when you pick clubs, diamonds or spades). The odds is simply the ratio of these two numbers, \(p/(1-p)\). In general: \[odds = \frac{p}{1-p}\] and \[p = \frac{odds}{odds+1}\] So the odds of picking hearts are \(0.25/0.75 \approx 0.333\), or \(1:3\). If you bet on hearts the odds are 1 to 3 in your disfavour. When two possible outcomes have the same probability (i.e. \(p = 0.5\)) the odds are \(1\). When \(p = 0\) the odds are also \(0\), and when \(p = 1\) the odds are \(\infty\).

The probability of picking a pip card is \(40/52 = 10/13 \approx 0.769\), and the odds are \((10/13)/(3/13) = 10/3 \approx 3.333\) or 10 to 3 in your favour. If you know the odds are 3.333 then \(p = 3.333/4.333 \approx 0.769\).

How odds of 1 to 3 is easier to grasp (for some) than a probability of 1 in 4, is one of the great mysteries of the world. Some believe it’s evidence for two different human species coexisting, possibly an extraterrestrial alien.

Statisticians sometimes use odds for mathematical reasons, and everybody hates them for it. But they have good reasons, as we shall see.

One of the situations where odds are useful is when modelling binary (1/0) either/or outcomes, such as success/failure or yes/no. We want to understand how the probability of, say, success changes with one or more explanatory variables.

Say we want to study if surviving the sinking of the Titanic was related to the passenger’s age. A normal simple linear regression would not be so good: we only have 1/0 response data, the errors cannot be normally distributed etc, and the regression may predict impossible values of larger than 1 or smaller than zero probabilities.

We could calculate the survival rate for different age groups and use this as our dependent variable instead of 1/0, but this is not ideal: the variance is not equal and not normally distributed, we lose the information about the sample size behind each proportion, and we run into problems if we want to include more explanatory variables. We need a better way to convert this into a continuous variable so that we can make use of the machinery of regression models. But thinking about the outcome as proportions or probabilities is a good start—we want to model the probability that survived = 1.

In a normal regression we model the expected value of the dependent variable as a linear function of the explanatory variable(s). What we can do here is to model the expected value of a function of the expected value of the dependent variable. This link function needs to respect that probabilities are bounded by \([0,1]\) while at the same time take on any value itself. This is where odds come in.

# Create a function to convert probabilities to odds:
probs.to.odds <- function (p) {
  stopifnot(p >= 0 & p <= 1)
  p/(1-p)}

# Check what the odds are for some probabilities:
probs.to.odds(c(0, 0.05, 0.5, 0.95, 0.99, 1))

[1]  0.00000000  0.05263158  1.00000000 19.00000000 99.00000000         Inf

# Make a sequence of probabilities from 0.001 to 0.999, incremented in steps of 0.001:
probabilities <- seq(0.001, 0.999, 0.001)

# Convert the probabilities to odds, using our function:
odds <- probs.to.odds(probabilities)

plot(x = probabilities, y = odds, ylim = c(0, 1)) # limit to odds under 1

plot(x = probabilities, y = odds, ylim = c(0, 100)) # limit to odds under 100

What if we take the logarithm of the odds?

Bingo!

plot(x = probabilities, y = log(odds)) 

# instead of log(odds) could have used the built-in function logit: 
# plot(y = probabilities, x = logit(probabilities)) 

Probability goes from 0 to 1, odds go from 0 to infinity, and by taking the logarithm of the odds we can get from \(-\infty\) to \(+\infty\) in a beautiful, symmetrical way.

\(\log(p/1-p)\) is so useful that it has not only one but two nicknames: the log-odds of an event is the logit of the probability of the event.

# check that the logit values are the same as the log(odds):
all.equal(logit(probabilities), log(odds))

[1] TRUE

The point of this was partly to escape the boundary of [0,1] for probabilities (since the log-odds can take on any real number), but there are other mathematical reasons too for using log-odds that has to do with ease of computation. There are also more nuances that can be taken care of by using a generalized linear model (glm). When we use the logit as the link function and binomial error distribution the glm model is called logistic regression.

In a logistic regression the coefficients are log(odds ratio). Yes, odds ratios, not odds. (Except for the intercept, which is the log(odds) when all other coefficients (if any) are at their reference levels — if no other coefficients this would simply be the overall log(odds) of success).

The coefficients are odds ratios (OR) because they are used to compare the odds between groups, we can express this as a ratio of two odds, i.e. a ratio of ratios. For example, if the probability of great success using EX is \(0.49\) and when not using EX is \(0.46\), the odds ratio is:

\[\frac{0.49/(1-0.49)}{0.46/(1-0.46)} = \frac{0.96078}{0.85185} = 1.128\] That is, using EX has 12.8% better odds of success, or 1.128 times better odds. (This does not mean that OR is the only possible, or necessarily the best, effect size measure for comparing proportions. This is a complex issue.)

So, to interpret logistic regression results we have to understand odds ratios. The actual coefficients are log(OR), but its customary to exponentiate to get back to the odds scale and sane people would also show results on the original (probability) scale, for example as in the prediction plots here (Fig. 11).

If we take the inverse logit of log-odds we get back to probabilities again. And the inverse of the logit is the well-known s-shaped logistic growth function! \(f(x) = 1/(1+e^{-x})\)

#the logistic function:
inverselogit <- function (x) {1/(1+exp(-x))} 

# check that the logistic function really gets us back to probability:
all.equal(inverselogit(log(odds)), probabilities) 

[1] TRUE

As an example of logistic regression, we can model the survival probability of passengers on the RMS Titanic that sank in 1912, as a function of sex, age and travel class (and first-order interactions). In total approximately 1500 of the 2,224 crew and passengers did not survive. We have complete data on these variables for 1046 of the 1316 passengers.

We skip tables of coefficients and only show the prediction plots here, on the probability scale (1 = survived, 0 = not survived). Included in the background of the plots are the age distributions of survivors and non-survivors, for males and females (as dot plots).

First a linear model:

lm(survived ~ sex+age+pclass+sex:age+age:pclass+sex:pclass, data = titanic)

Then a logistic regression model:

glm(survived ~ sex+age+pclass+sex:age+age:pclass+sex:pclass, family = binomial(link = "logit"), data = titanic)

We see that this kind of binary outcome data is handled poorly by the normal linear regression model, making nonsensical predictions of more than 1 and less than 0 probability of survival (although it does capture some of the same features as the glm).

Only four female first-class passengers died on the Titanic, one because she refused to leave her Great Dane. Three of the 12 dogs on board survived. They were all travelling first class.

Success rate: The EX factor

FIGURE 9: Proportion great succes in the different treatment groups, with binomial confidence intervals and sample sizes.

On land (a)

At sea (b)

Combined (c)

FIGURE 10: Proportion great success in treatments with and without special production EX. Production done at sea was excluded in a) because workshop ship skill (only effective at sea) vs. Meister title with nearby apprentice (only effective on land) may differ slightly in the respective boost to success rate they give. b) shows the results at sea, and c) the combined results.

The rate of great success did not vary dramatically among the treatment groups (Fig. 9). Looking separately at production done at land and sea (but ignoring refining and skill rank), in both cases crafting with the EX equipped gave slightly better success than without EX (+3.1 % on land and +4.4% at sea, Fig. 10). Similarly, production at sea gave better success than producing on land, again not by much.

This latter comparison is in fact an estimate of the difference in success rate from being in ship with workshop at sea vs being on land with Meister title and apprentice nearby, and has probably nothing to do with the sailor equipment. These two effects are mutually exclusive since the Meister effect is only active on land and workshop only at sea. The data (compare Fig. 10 a and b) suggest a difference in success rate between land and sea of +2.1% (without EX) or +3.4% (with EX). According to http://gvo.gamedb.info/wiki/ the effect of Meister title with apprentice is +5% and of workshop +10%, so I had expected the difference to have been slightly larger, +5%. However, 5% is still well within the uncertainty of the estimates.

TABLE 3: Logistic regression model of the effect of special production EX and producing at sea (with workshop) on success rate. See Tab. 6 in the Appendix for other models considered.

Variable	Great Success Rate	OR1	95% CI	p-value
EX				0.012
0	1,223 / 2,600 (47.0%)	—	—
r6	1,313 / 2,600 (50.5%)	1.15	1.03–1.28
where				0.067
land	1,745 / 3,640 (47.9%)	—	—
sea	791 / 1,560 (50.7%)	1.12	0.99–1.26
Abbreviations: CI = Confidence Interval, OR = Odds Ratio
1 Contrasts: dummy-coding

Statistical modelling confirmed what looking at the raw data told us (see Tab. 3 and Tab. 6), and suggests that the effect of the EX is similar, or slightly larger, than the effect of land/sea (with considerable overlap in the uncertainty). The odds ratios were 1.15 and 1.12 respectively, which means that the odds of a great success is 15% better when producing with EX on compared to without, and 12% better when producing at sea with workshop compared to at land with meister title. These effects are additive (i.e. they stack), so on land without EX there was about 46% probability of great success and at sea with EX about 52% (Fig. 11).⁸

⁸  Of course, the other factors that were held constant such as Oxford production success 1-4, paymaster 120SS aide and The Magician tarot card would also affect the actual probabilities.

FIGURE 11: Predicted probabilities of great success from the logistic regression model in Tab. 3, with 95% confidence intervals.

Degreee of success: The SE factor

Piercing Power

Strip plots (a)

Barplots (b)

Interaction plots (c)

Alternative ways to display the factors

FIGURE 12: Piercing power of minion(4) cannons made in the different treatment groups. a) Individual data poins, scaled by the number of observations. Also shown is the mean values and their confidence intervals obtained by bootstrap resampling, as well as the overall mean across all observations. The values on top represent the sample size in each group. b) shows the same data but as barplots (frequency of values within each treatment group). c) shows the mean values in the groups, ignoring EX.

FIGURE 13: Predicted piercing power from the model in Tab. 4 with 95% confidence intervals.

Piercing power did not vary much among the treatment groups (Fig. 12 a). The distribution of PP was somewhat skewed (Fig. 12 b); approximately 1/3 of the masters’s minions produced were tier 1 and 1/6 in each of the remaining 4 tiers, irrespective of the level of the manipulated factors. It is possible that there was a very small (around 0.5 PP, or 1/8 PP tier) effect of sailor equipment (Fig. 13, but the evidence was not strong (Tab. 4). None of the other factors explained any variation at all (see Tab. 7 in the Appendix).

TABLE 4: Linear model of the effect of sailor equipment on piercing power (with only one factor, and only two levels, this is equivalent to a t-test). The evidence is only weak, for a very small effect. See Tab. 7 in the Appendix for other models considered, and Fig. 13.

Effect	df	F	pes1	p
sailor equip	1, 1546	2.55	0.002	0.110
1 partial eta-squared (effect size).

Durability

The results for durability were spectacularly different than for piercing power (Fig. 14).

Durability of any master’s cannon range from 105% to 140% of the base (100 base for master’s minion 4) and, unlike piercing power, can take on any integer value in this range. Having refined the skill increased the minimum to 110%, as did the sailor equipment (Fig. 14). However, the increased minimum improvement far from explain all that is going on in Fig. 14.

Strip plots (a)

Histograms (b)

Interaction plots (c)

Alternative ways to display the factors

FIGURE 14: Durability of minion(4) cannons made in the different treatment groups. a) Individual data points, scaled by the number of observations. Also shown is the mean values and their confidence intervals obtained by bootstrap resampling, as well as the sample size in each group. b) shows the same data, but as histograms (binwidth = 8). c) shows the mean values in each group, ignoring EX.

TABLE 5: Anova table (Type III tests) for the chosen linear model of effects on durability (full factorial model but no effect of EX). See Tab. 8 in the Appendix for other models considered, and Fig. 15 for model predictions.

Effect1	df	F	pes2	p
sailor equip	1, 1540	296.78	0.162	<0.001
refined	1, 1540	257.20	0.143	<0.001
crafting rank	1, 1540	151.10	0.089	<0.001
sailor equip x refined	1, 1540	128.97	0.077	<0.001
sailor equip x crafting rank	1, 1540	16.57	0.011	<0.001
refined x crafting rank	1, 1540	29.36	0.019	<0.001
sailor equip x refined x crafting rank	1, 1540	4.13	0.003	0.042
1 x indicates interaction. 2 partial eta-squared (effect size).

The first thing to notice is how much lower the values are with unrefined compared to refined skill. Second, producing with the sailor equipment on gives much larger values than without, especially for unrefined and low skill rank. Third, increasing the skill rank from r3 to r20 seems to have a large impact, especially with unrefined skill. Furthermore, there does not seem to be any effect of the EX (compare the four pairs that differ only in EX in Fig. 14).

The above was confirmed by statistical modelling (Tab. 5, Fig. 15). There was no effect of the EX and it was dropped from the model, but all other factors and their interactions had to be kept and together account for 44% of the variation in durability (see Tab. 8 in the Appendix for how this model was selected).

All the factors interact with each other so it is not meaningful to give overall estimates of how large the effect is of either skill rank, refining or using SE: the effect depends on the level of the other factors.

Even at maximum rank and having refined the skill, using the sailor equipment still yields better results than without (Fig. 15). However, the effect of the SE is much larger with unrefined skill — r20 unrefined with SE was in fact indistinguishable from refined r20 with SE (or only very slightly lower), although I suspect there may be a larger difference for a higher rank recipe.

There was no correlation between PP and durability (r_s = −0.029); these parameters are set independently from each other and there is no overall tier of success.

FIGURE 15: Predicted durability from the linear model in Tab. 5, with 95% confidence intervals.

Conclusions: are the special production SE and EX any good?

The SE and EX stack, in the sense that they are not for the same effect: It seems clear that the EX increases the rate of success and the SE the degree of success, as the descriptions indicated (Fig. 1, Fig. 2). But this comes with some important caveats, as discussed below.

The SE

It turned out to be important to include durability in this investigation, and test with unrefined and low rank casting skill. This made it obvious there was a huge effect of the sailor equipment (Fig. 15), but mostly at low/unrefined rank and only on durability, not piercing power (Fig. 13). At the same time, since durability is of little concern and most people would produce at maximum refined rank, this revealed the effect to be rather trivial in the case of cannons. So, yeah, great success with the SE is more likely to be “high-performing” as the description of the SE advertises (Fig. 2) but not for what matters (piercing power).

Usually in UWO, effects are not a relative percentage but instead comes in the form of a fixed added number of percentage points. Such as +30% grading success. This is definitely not the case for degree of success from refining, skill rank or sailor equipment. These effects are not additive at all.

The lack of any effect on piercing power is likely not due to the low rank of the Minion recipe, as there was no effect of skill rank or refining either. However, this should be directly tested on high rank cannons to be sure.

As revealed here, forging is not considered “special production” and we have previously shown that the SE has no effect in forging recipes⁹. However, forging has its own exclusive skill with a large effect (item improvement techniques from Oxford) so one could say that the SE is “Oxford for non-forging” but unfortunately, at least for cannons, the effect is so small (or absent) for PP that I think it misses its mark.

⁹ Stinker & Marsali (2024) Forging success in UWO — experiments with fur boots

¹⁰ Skill rank was not tested with the Gorz.

¹¹ Stinker (2024) Crafting success in UWO — experiments with Gorz

Could the SE have more relevant effects on other equipment than cannons? It might be that the SE only affects the durability, and not attack or defence, on personal gear as well. However, whereas refining and skill rank had no effect on PP on the cannons (only on durability), refining¹⁰ had a large effect on both attack, defence and durability on the Gorz weapon¹¹. I suspect therefore that the SE will indeed increase attack/defence of personal gear and not only durability. Future studies should investigate this.

The EX

The EX¹² had no effect on the piercing power or durability of master’s minions, but a small effect on the rate of great success of around +3% (Fig. 11). Since the effect it so small it required very large sample sizes to pick up a reliable signal in all the random noise.

¹² The golden arc is not the only special goods production EX r6. There is also ring of royalty and custom amber ring, as well as the r4 ancient priest’s knowledge book.

¹³ Stinker (2024) Crafting success in UWO — experiments with Gorz

It seems I may have dismissed the EX effect in the earlier Gorz study too quickly¹³. This was partly due to one of the comparisons in that study happened to give slightly better success without than with the EX, and a small effect overall, not “significantly” different from zero. But in fact, the data on Gorz is as compatible with an effect of +5% as with an effect of 0% and the observed difference in great success rate is not far from that seen with minions: 47.8% with EX and 45.6% without EX = 2.2%. For minions (on land): 46.4% vs 49.5 = 3.1% and at sea: 48.5% vs 52.9% = 4.4%

It is also of course possible that fixed recipes at NPCs such as the Gorz does not count as special goods production (analogous to how forging does not count), but this would be very hard to confirm given the small effect of the EX and the impossibility to test such recipes with the sailor equipment. However, I see no reason why such recipes should not count and the most parsimonious explanation given all the currently available evidence is that the EX did work with the Gorz in the same way as with the cannons. It may also suggest that the low effect found in the present study is is not due to the low required skill rank since Gorz is a r18 recipe.

No effect of EX was found in alchemy experiments¹⁴, but this is as expected since the SE explicitly states that it is restricted to sewing, casting and handicrafts. This is not stated for the EX but it is reasonable to assume there is only one definition of “special goods production”.¹⁵ The EX effect was also absent from the study of forging¹⁶, but as revealed here forging is not considered to be “special goods production” either.

¹⁴ Stinker (2023) Report from a FABtory

¹⁵ Although, there is also “Special Production” Oxford research which confusingly enough refers to production with alchemy/linguistics, and “Special Prod. Tech” which is Oxford research that refers to Florence production. Makes me wonder, how special can you get?

¹⁶ Stinker & Marsali (2024) Forging success in UWO — experiments with fur boots

We are therefore limited to casting, sewing and handicrafts recipes that are not forging and effectively not cannons, and possible to do at sea for the SE, and not worth using master’s secret on (for the EX). It is not that many recipes left.

---

In summary, the EX effect is small, and this is why it was missed earlier in a study of crafting weapons. The SE is perhaps more interesting than the EX in that it can have huge effects. But only at low crafting rank and/or unrefined crafting skill which in practice is not a concern. Also disappointing is that the SE effect was only for the durability and not the piercing power of cannons. Higher rank recipes and recipes different from cannons should be tested in future studies. The potential for the SE seems somewhat limited, as most of the best recipes are fixed at NPCs on land and thus the sailor equipment cannot be used for these. Further, great success rate (and thus the EX) can be made irrelevant by using master’s secrets.

Acknowledgements

Thanks to Marsali and TardyMcGormless for discussions.

Appendix

Statistical Models

Model selection, coefficients and model predictions

Rate of success

Model selection

TABLE 6: Logistic regression models for the rate of great success. Model 6 was selected (additive effect of EX and location). Model 5 had somewhat lower AIC but is rejected, because the odds ratio suggesting possible lower succcess for r20 than r3 casting is a priori very unlikely. All production at sea were done with sailor equipment, but note that location (“where”) is here (for success rate) interpreted as the effect of workshop at sea vs meister title on land, and not as an effect of sailor equipment per se. The models are ranked according to AIC in a separate table under the main table.

Parameters2	Candidate Models (logistic regression)1
	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
(Intercept)	0.88 [0.69, 1.13] p = 0.321	0.94 [0.79, 1.12] p = 0.483	0.94 [0.81, 1.09] p = 0.407	0.87 [0.79, 0.95] p = 0.002	0.93 [0.82, 1.05] p = 0.226	0.86 [0.79, 0.94] p = <0.001	0.92 [0.86, 0.98] p = 0.013	0.89 [0.82, 0.96] p = 0.003	0.95 [0.90, 1.01] p = 0.076
EX [0]	—	—	—	—	—	—		—
EX [r6]	1.43 [1.01, 2.02] p = 0.044	1.34 [0.99, 1.81] p = 0.055	1.16 [1.03, 1.32] p = 0.014	1.13 [0.99, 1.29] p = 0.068	1.17 [1.05, 1.31] p = 0.005	1.15 [1.03, 1.28] p = 0.013		1.15 [1.03, 1.28] p = 0.013
where [land]	—	—	—	—	—	—	—
where [sea]	1.18 [0.84, 1.67] p = 0.335	1.06 [0.83, 1.35] p = 0.664	1.08 [0.95, 1.22] p = 0.240	1.08 [0.92, 1.28] p = 0.341	1.08 [0.96, 1.23] p = 0.213	1.12 [0.99, 1.26] p = 0.067	1.12 [0.99, 1.26] p = 0.068
crafting_rank [r3]	—	—	—		—
crafting_rank [r20]	0.84 [0.64, 1.10] p = 0.196	0.89 [0.73, 1.09] p = 0.271	0.91 [0.78, 1.06] p = 0.221		0.89 [0.79, 1.02] p = 0.084
refined [no]	—		—
refined [yes]	1.13 [0.80, 1.60] p = 0.483		0.97 [0.83, 1.13] p = 0.717
EX [r6] × where [sea]	0.68 [0.42, 1.11] p = 0.121	0.76 [0.50, 1.16] p = 0.209		1.06 [0.84, 1.34] p = 0.629
EX [r6] × crafting_rank [r20]	0.81 [0.52, 1.25] p = 0.345	0.83 [0.60, 1.16] p = 0.286
where [sea] × crafting_rank [r20]	1.07 [0.69, 1.65] p = 0.760	0.95 [0.66, 1.37] p = 0.796
EX [r6] × refined [yes]	0.98 [0.64, 1.52] p = 0.945
where [sea] × refined [yes]	0.79 [0.49, 1.29] p = 0.352
EX [r6] × where [sea] × crafting_rank [r20]	1.76 [0.92, 3.39] p = 0.090	1.65 [0.96, 2.82] p = 0.070
EX [r6] × where [sea] × refined [yes]	0.94 [0.49, 1.81] p = 0.850
Num.Obs.	5200	5200	5200	5200	5200	5200	5200	5200	5200
AIC	7208.9	7203.5	7202.9	7203.8	7201.0	7202.0	7206.2	7203.3	7207.6
BIC	7287.6	7256.0	7235.7	7230.0	7227.2	7221.7	7219.4	7216.5	7214.1
1 Values shown are the Odds Ratios and [95% confidence intervals]. Contrasts: dummy coding.
2 × indicates interaction

Model	K	AICc	ΔAICc
(5)	4	7201.03	0.00
(6)	3	7202.00	0.98
(3)	5	7202.90	1.87
(8)	2	7203.35	2.32
(2)	8	7203.57	2.54
(4)	4	7203.77	2.75
(7)	2	7206.24	5.21
(9)	1	7207.58	6.55
(1)	12	7208.99	7.96

Model predictions

FIGURE 16: Predicted probabilities of great success from model (6), with 95% confidence intervals.

Piercing Power

Model selection

TABLE 7: Linear models for piercing power. Model 5 was selected (effect only of sailor equipment). The null model (model 6) is almost as good, which means there is very litte evidence of any effects at all. Sailor equipment explains almost nothing of the variation in piercing power. The models are ranked according to AIC in a separate table under the main table.

Parameters2	Candidate Models1
	(1)	(2)	(3)	(4)	(5)	(6)
(Intercept)	158.68 [158.36, 158.99] p = <0.001	158.66 [158.36, 158.96] p = <0.001	158.66 [158.36, 158.96] p = <0.001	158.66 [158.36, 158.96] p = <0.001	158.66 [158.36, 158.95] p = <0.001	158.66 [158.37, 158.96] p = <0.001
Sailor Equip1	-0.31 [-0.62, 0.01] p = 0.058	-0.24 [-0.54, 0.06] p = 0.114	-0.24 [-0.54, 0.06] p = 0.113	-0.24 [-0.54, 0.05] p = 0.108	-0.24 [-0.54, 0.05] p = 0.110
Refined1	0.02 [-0.34, 0.39] p = 0.906	0.03 [-0.33, 0.39] p = 0.869
Crafting Rank1	-0.12 [-0.49, 0.24] p = 0.504	-0.13 [-0.50, 0.23] p = 0.474	-0.12 [-0.43, 0.20] p = 0.462
Ex Rank1	0.14 [-0.23, 0.50] p = 0.467	0.13 [-0.23, 0.50] p = 0.468	0.12 [-0.20, 0.43] p = 0.460	0.08 [-0.22, 0.38] p = 0.599
Sailor Equip1 × refined1	0.02 [-0.30, 0.34] p = 0.909
Sailor Equip1 × crafting Rank1	0.02 [-0.30, 0.34] p = 0.898
Refined1 × crafting Rank1	-0.06 [-0.38, 0.26] p = 0.708
Sailor Equip1 × refined1 × crafting Rank1	0.21 [-0.11, 0.52] p = 0.200
Num.Obs.	1548	1548	1548	1548	1548	1548
R2	0.003	0.002	0.002	0.002	0.002	0.000
R2 Adj.	-0.002	-0.000	0.000	0.001	0.001	0.000
AIC	9925.2	9919.1	9917.2	9915.7	9914.0	9914.6
BIC	9978.7	9951.2	9943.9	9937.1	9930.0	9925.2
Std.Errors	HC3	HC3	HC3	HC3	HC3	HC3
1 Contrasts: effects coding. Variance-covariance matrix: HC3 (standard errors robust to heteroskedasticity)
2 × indicates interaction

Model	K	AICc	ΔAICc
(5)	3	9914.01	0.00
(6)	2	9914.56	0.55
(4)	4	9915.75	1.73
(3)	5	9917.22	3.20
(2)	6	9919.20	5.19
(1)	10	9925.38	11.37

Model predictions

FIGURE 17: Predicted piercing power from model (5), with 95% confidence intervals.

Durability

Model selection

TABLE 8: Linear models for durability. Model (3) was selected (full factorial model but no effect of EX). It is possible that the higher order interaction could also be be removed (model 5, ΔAICc = 2.12), but the second-order interactions and the main effects (apart from EX) must clearly be kept. The evidence is very strong for large effects of these factors, together they account for 44% of the total variation in durability. But the presence of interactions means we cannot meaningfully say how large the main effects are because they depend on the level of other factors, and interpretation is best done from plots (Fig. 15). The models are ranked according to AIC in a separate table under the main table.

Parameters2	Candidate Models1
	(1)	(2)	(3)	(4)	(5)	(6)
(Intercept)	132.30 [131.84, 132.76] p = <0.001	131.88 [131.42, 132.34] p = <0.001	132.30 [131.84, 132.76] p = <0.001	131.87 [131.41, 132.34] p = <0.001	132.29 [131.83, 132.75] p = <0.001	131.80 [131.23, 132.38] p = <0.001
Sailor Equip1	-4.02 [-4.49, -3.56] p = <0.001	-4.29 [-4.76, -3.82] p = <0.001	-4.02 [-4.49, -3.56] p = <0.001	-4.29 [-4.76, -3.82] p = <0.001	-4.18 [-4.61, -3.74] p = <0.001
Refined1	-3.71 [-4.25, -3.18] p = <0.001	-3.58 [-4.15, -3.02] p = <0.001	-3.74 [-4.21, -3.28] p = <0.001	-3.71 [-4.19, -3.23] p = <0.001	-3.77 [-4.23, -3.30] p = <0.001
Crafting Rank1	-2.90 [-3.44, -2.36] p = <0.001	-3.02 [-3.58, -2.46] p = <0.001	-2.87 [-3.33, -2.41] p = <0.001	-2.90 [-3.38, -2.42] p = <0.001	-2.85 [-3.32, -2.39] p = <0.001
Ex Rank1	0.06 [-0.48, 0.60] p = 0.824	0.25 [-0.34, 0.83] p = 0.405
Sailor Equip1 × refined1	-2.65 [-3.11, -2.19] p = <0.001		-2.65 [-3.12, -2.19] p = <0.001		-2.65 [-3.12, -2.19] p = <0.001
Sailor Equip1 × crafting Rank1	-0.95 [-1.41, -0.48] p = <0.001		-0.95 [-1.41, -0.49] p = <0.001		-0.96 [-1.43, -0.49] p = <0.001
Refined1 × crafting Rank1	-1.26 [-1.73, -0.80] p = <0.001		-1.26 [-1.73, -0.80] p = <0.001		-1.25 [-1.71, -0.79] p = <0.001
Sailor Equip1 × refined1 × crafting Rank1	-0.47 [-0.94, -0.01] p = 0.045		-0.47 [-0.94, -0.01] p = 0.045
Num.Obs.	1548	1548	1548	1548	1548	1548
R2	0.443	0.358	0.443	0.358	0.441	0.000
R2 Adj.	0.440	0.357	0.440	0.357	0.439	0.000
AIC	11092.2	11302.6	11090.2	11301.3	11092.4	11981.3
BIC	11145.6	11334.7	11138.3	11328.1	11135.1	11992.0
Std.Errors	HC3	HC3	HC3	HC3	HC3	HC3
1 Contrasts: effects coding. Variance-covariance matrix: HC3 (standard errors robust to heteroskedasticity)
2 × indicates interaction

Model	K	AICc	ΔAICc
(3)	9	11090.33	0.00
(1)	10	11092.31	1.98
(5)	8	11092.46	2.12
(4)	5	11301.37	211.04
(2)	6	11302.66	212.33
(6)	2	11981.31	890.98

Model predictions

FIGURE 18: Predicted durability from model (3), with 95% confidence intervals.

Batch variation

Note 2: Permutation tests

The most important assumption in statistical analyses of data is commonly that the observations are independent and identically distributed, that is that they are a random sample. Biases in the random number generator used in UWO is often claimed, and would violate this fundamental assumption.

When crafting it is often one gets the feeling that success or certain values are clustered. Some batches may yield very many great success, or mostly (or only!) maximum durability or minimum PP, for example.

To test if the variation within batches differ from what can be expected from random variation, I calculated the standard deviation (SD) of piercing power and durability for each batch of 20 cannons, and compared the mean SD of the 156 batches with what can be expected from random variation using permutation tests.

The random variation was simulated by reshuffling (sampling without replacement) the values of the 1548 great success cannons among batches and then recalculating the SD values on the permuted sample. Since treatments could influence the values, this was done stratified (i.e. the permutations were done within each treatment). This procedure was then repeated 10,000 times to get an expected distribution of SD values, and the probability of obtaining the observed value could then be obtained from this distribution.

Normal success was excluded when calculating both the observed and simulated SD values. This meant that batches would have less than 20 cannons, and the simulations therefore used the same numbers in each batch as the observed.

This exercise showed that the variation between batches was entirely as expected from random variation, both in piercing power and in durability (figures below).

I have, over this series of investigations, used a suite of different statistical techniques to verify that the RNG algorithms used when crafting in UWO are unbiased and indistinguishable from a random process. This has included time-series analyses of auto-correlation; runs-tests; randomisation tests; and now permutation tests.

Next I want to confirm that grading success is entirely as expected as well. If anyone has collected data (or want to) on this grandmother of all painful data gathering and want to collaborate, let me know.

R code used in this report

R code

knitr::opts_knit$set(global.par = TRUE)
library(ggplot2)
library(dplyr)
library(ggh4x)
library(legendry)
library(ggtext)
library(gt)
library(gtsummary)
library(afex)
library(tidyr)
library(purrr)
library(boot)
library(ggstats)
library(cowplot)
library(modelsummary)
library(sandwich)
library(ggdist)
library(patchwork)

options(pillar.sigfig = 5)
par(pty="m", mar =c(3.7, 3.5, 0.5, 0.5), mgp= c(2.5, 1, 0))
config_modelsummary(factory_default = "gt")

minions <- as.data.frame(readODS::read_ods("minion4.ods", sheet =1))
minions$treatment <- factor(minions$treatment)
minions$land_sea <- factor(minions$land_sea)
minions$workshop <- factor(minions$workshop)
minions$sailor_equip <- factor(minions$sailor_equip, levels = c("(land)", "r5"))
minions$sailor_equip[is.na(minions$sailor_equip)] <- "(land)"
minions$refined <- factor(minions$refined)
minions$crafting_rank <- factor(minions$crafting_rank, levels = c("r3", "r20"))
minions$Ex_rank <- factor(minions$Ex_rank)
minions$outcome <- factor(minions$outcome, levels = c("normal", "gs"))
minions$placeholder <- c(NA)
minions$placeholder <-  factor(minions$placeholder, levels = c("NA"))
summary(minions)

overalProp <- minions |>
  summarise(
    normal = sum(outcome == "normal"),
    gs = sum(outcome == "gs"),
    total = n()
  ) |>
  mutate(
    propGs = gs / total,
    propNormal = normal / total,
  )

minionsProp <- minions |>
  group_by(treatment) |>
  summarise(normal = sum(outcome == "normal"),
            gs = sum(outcome == "gs"),
            total = n(),
            land_sea = unique(land_sea),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank),
            oxford_item = unique(oxford_item),
            Ex_rank = unique(Ex_rank),
            placeholder = unique(placeholder))

additional <- as.data.frame(readODS::read_ods("minion4.ods", sheet =2))
additional$land_sea[is.na(additional$land_sea)] <- "land"
additional$land_sea <- factor(additional$land_sea, levels = c("land", "sea"))
additional$sailor_equip <- factor(additional$sailor_equip, levels = c("(land)", "r5"))
additional$sailor_equip[is.na(additional$sailor_equip)] <- "(land)"
additional$placeholder <- c(NA)

addidtionalProp <- additional |>
    group_by(treatment) |>
  summarise(normal = sum(nNormal),
            gs = sum(nGs),
            total = sum(sum),
            land_sea = unique(land_sea),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank),
            oxford_item = unique(oxford_item),
            Ex_rank = unique(Ex_rank),
            placeholder = unique(placeholder))

combinedProp <- rbind(minionsProp, addidtionalProp) |>
      group_by(treatment) |>
  summarise(normal = sum(normal),
            gs = sum(gs),
            total = sum(total),
            land_sea = unique(land_sea),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank),
            oxford_item = unique(oxford_item),
            Ex_rank = unique(Ex_rank),
            placeholder = unique(placeholder))

#binomial confidence intervals
binomCIgs <- binom::binom.wilson(combinedProp$gs, combinedProp$total)

CIgs <- cbind(combinedProp, binomCIgs)

#restructure for outcome and numbers, need treatment and #factoreroutcome (normal gs) and n
pivoted <- CIgs |>
  dplyr::select(treatment:placeholder) |>
  pivot_longer(cols = c(normal, gs), names_to = "outcome", values_to = "antall")
pivoted$outcome <- factor(pivoted$outcome, levels = c("normal", "gs"))

##slå sammen de med ex og de uten, men kun på land! og lag nye CIs
EXcombinedProp <- combinedProp |>
  filter(land_sea == "land") |>
  group_by(Ex_rank) |>
  summarise(normal = sum(normal),
            gs = sum(gs),
            total = sum(total))
EXbinomCIgs <- binom::binom.wilson(EXcombinedProp$gs, EXcombinedProp$total)
EXCIgs <- cbind(EXcombinedProp, EXbinomCIgs)
EX_X2 <- prop.test(EXcombinedProp$gs, EXcombinedProp$total, alternative = "less")

EXpivoted <- EXCIgs |>
  dplyr::select(Ex_rank:total) |>
  pivot_longer(cols = c(normal, gs), names_to = "outcome", values_to = "antall")
EXpivoted$outcome <- factor(EXpivoted$outcome, levels = c("normal", "gs"))

meanGsEX <- minions |>
  group_by(treatment) |>
  summarise(mean_pp = mean(pp, na.rm = TRUE),
            sd_pp = sd(pp, na.rm = TRUE),
            mean_dura = mean(durability, na.rm = TRUE),
            sd_dura = sd(durability, na.rm = TRUE),
            n = n(),
            valid_n = sum(!is.na(pp)))

# bootstrapped confidence intervals for pp and defence
### function needed for the bootstrap:
mean.fun <- function(x, index){mean(x[index], na.rm = TRUE)}

obsMeansPp <- minions |>
  group_by(treatment) |>
  summarise(mean_pp = mean(pp, na.rm = TRUE),
            sd_pp = sd(pp, na.rm = TRUE),
            n_pp = sum(!is.na(pp)),
            sailor_equip = unique(sailor_equip),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank),
            oxford_item = unique(oxford_item),
            Ex_rank = unique(Ex_rank),
            placeholder = unique(placeholder))

treatments <- unique(obsMeansPp$treatment)
n_treatments <- length(unique(obsMeansPp$treatment))
### set up dataframe to collect results for plotting:
pldatPp <-cbind(obsMeansPp, mean=rep(NA,n_treatments),
                low=rep(NA, n_treatments),
                upp=rep(NA, n_treatments))
rownames(pldatPp) <-c(treatments)
pldatPp

set.seed(666)
for(i in c(treatments)){
### the boot sample and the ci's:
    boot<-boot(minions$pp[minions$treatment==i], mean.fun, R=10000)
    ci<-boot.ci(boot,type = c("norm", "basic", "perc", 
    "bca"))
### get ci's (method: bca)
    ci[2]->pldatPp[i,"mean"] # mean
    ci$bca[1,4]->pldatPp[i,"low"]  # lower ci
    ci$bca[1,5]->pldatPp[i,"upp"]} # upper ci
pldatPp

set.seed(666)
obsMeansDur <- minions |>
  group_by(treatment) |>
  summarise(mean_dur = mean(durability, na.rm = TRUE),
            sd_dur = sd(durability, na.rm = TRUE),
            n_dur =  sum(!is.na(durability)),
            sailor_equip = unique(sailor_equip),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank),
            oxford_item = unique(oxford_item),
            Ex_rank = unique(Ex_rank),
            placeholder = unique(placeholder))

pldatDur <-cbind(obsMeansDur, mean=rep(NA,n_treatments),
                 low=rep(NA, n_treatments),
                 upp=rep(NA, n_treatments))

rownames(pldatDur) <-c(treatments)
pldatDur

for(i in c(treatments)){
    boot<-boot(minions$durability[minions$treatment == i], mean.fun, R = 10000)
    ci<-boot.ci(boot,type = c("norm", "basic", "perc", 
    "bca"))
    ci[2] -> pldatDur[i, "mean"] # mean
    ci$bca[1, 4]->pldatDur[i, "low"]  # lower ci
    ci$bca[1, 5]->pldatDur[i, "upp"]
    } # upper ci
pldatDur

#to check what the index values are for upper and lower ci in boots.ci. Can delete
# test <- minions |>
#   filter(treatment == "D")
# test2 <- boot(test$durability, mean.fun, R=1000)
# test3 <- boot.ci(boot, type = c("norm", "basic", "perc", 
#     "bca"))
# test3
# test2

png <- magick::image_read("images\\minion4.png")

img <- grid::rasterGrob(png, interpolate = TRUE)

duraLoss <- as.data.frame(readODS::read_ods("minion4.ods", sheet =4))
specialTable <- as.data.frame(readODS::read_ods("minion4.ods", sheet =3))

reducedMinions <-minions[!is.na(minions$pp), ]

#showtext::showtext_auto()


theme_set(theme_half_open(font_size = 12, font_family = "Roboto")) +
theme_update(plot.background = element_rect(fill = "grey92", 
                                            colour = NA), 
             legend.background = element_rect(fill = "grey92"),
             legend.key = element_rect(linewidth = 0.5),
             legend.key.spacing.y = unit(1.5, "pt"),
             legend.margin = margin(0, 0, 0, 0),
             plot.margin = margin(c(7, 7, 7, 7), unit = "pt"))
update_geom_defaults("text", list(family = theme_get()$text$family))
update_geom_defaults("richtext", list(family = theme_get()$text$family))

update_geom_defaults("bar", list(fill = I("#377eb8"), 
                                 colour =("grey30"),
                                 linewidth = 0.3))

windowsFonts(sans = windowsFont("Roboto"))

specialTable |> dplyr::select(special_production:dura_loss) |>
  gt(rowname_col = "recipe_type", groupname_col = "special_production", row_group_as_column = F) |> 
  tab_spanner(label = md("Recipe"), 
               columns = c(2:4)) |>
   tab_spanner(label = md("Product"), 
               columns = c(5:6)) |>
   tab_spanner(label = md("Sailor Equip."), 
               columns = c(7:8)) |>
  tab_stubhead(label = md("Type")) |>
  sub_missing() |>
  cols_label(recipe_type = "Recipe type",
             normal = md("Normal success"),
             great_success = md("Great success"),
             dura_message = md("Message"),
             dura_loss = md("Loss"),
             skill = md("Skill"),
             rank = md("Rank")) |>
  tab_style(style = cell_text(color = "#777"), locations = cells_title ()) |>
  cols_align(align = 'center', columns = 7:8) |> 
  opt_stylize(style = 3, color = "blue") |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
  tab_footnote(footnote = "Denotes if great success resulted in a differently named product than normal success, or same product with better stats, or in larger quantity",
                locations = cells_column_labels(columns = great_success)) |>
  tab_footnote(footnote = "Durability loss. This is per great success, and comes in addition to the normal loss over time at sea",
                 locations = cells_column_labels(columns = dura_loss)) |>
  tab_style(style = cell_text(size = px(15), weight = "bold"),
            locations = cells_row_groups()) |>
  tab_style(style = cell_text(size = px(12)),
            locations = cells_footnotes()) |>
  tab_options(data_row.padding = px(2),
              summary_row.padding = px(3),
              grand_summary_row.padding = px(3),
              row_group.padding = px(4),
              footnotes.padding = px(2),
              quarto.use_bootstrap = TRUE) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_body()) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_stub())  |>
  tab_stub_indent(rows = everything(), indent = 3) |>
  opt_table_font(font = list("roboto"))

duraLoss |>
  ggplot(aes(x = reqSkillRank, y = duraLoss)) +
  geom_line() +
  geom_point(aes(colour= cat, alpha = cat), size = 4) +
  scale_colour_manual(values = c("#377eb8", "grey"))+
  scale_x_continuous(breaks = seq(2, 20, 2), limits =c(1,20))+ 
  scale_y_continuous(n.breaks = 6, limits = c(1,10)) +
  ylab("Durability Loss") +
  xlab("Required Skill Rank") +
 scale_alpha_manual(name = "", values = c(1, 0)) +
  theme(legend.position = "none") 

treatmentGroups <- minions |>
  group_by(treatment) |>  
  summarise(
    treatment = factor(unique(treatment)),
    Where = factor(unique(land_sea)),
    SE = recode(factor(unique(sailor_equip)), "(land)" = NA_character_,
                "r5" = "r5"),
    Refined = unique(refined),
    Casting = factor(unique(crafting_rank)),
    EX = recode(unique(Ex_rank),
                "0" = "no",
                "6" = "r6"),
    n_normal = sum(outcome == "normal"),
    n_gs = sum(outcome == "gs")) 
treatmentGroups$Extra <- c(1040, 1040, rep(0, 10))
treatmentGroups$Total <- treatmentGroups$n_normal+treatmentGroups$n_gs+treatmentGroups$Extra


treatmentGroups |>
gt(rowname_col = "treatment") |>
  tab_spanner(label = md("Factors"), 
              columns = c(3:6)) |>
  tab_spanner(label = md("Outcome"), 
              columns = c(7:8)) |>
  tab_stubhead(label = "Treatment") |>
  sub_missing() |>
  cols_label(SE = "Sailor Equip",
             n_normal = "Normal",
             n_gs = "GS") |>
  #tab_style(style = cell_text(color = "#777"), locations = cells_title ()) |>
  opt_stylize(style = 3, color = "blue") |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
  grand_summary_rows(
    columns = c(7:10),
    fns = list(Total ~ sum(.))) |>
  tab_footnote(footnote = "GS = Great Success. This also gives the sample sizes for piercing power and durability.",
               locations = cells_column_labels(columns = n_gs)) |>
  tab_footnote(footnote = "Additional sample only used for success rate with and without EX",
               locations = cells_column_labels(columns = Extra)) |>
  tab_style(style = cell_text(size = px(12)),
            locations = cells_footnotes()) |>
  tab_options(data_row.padding = px(2),
              summary_row.padding = px(3),
              grand_summary_row.padding = px(3),
              row_group.padding = px(4),
              footnotes.padding = px(2),
              quarto.use_bootstrap = TRUE) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_body()) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_stub())  |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_grand_summary())  |>
    tab_style(style = cell_text(size = px(15)),
            locations = cells_stub_grand_summary())  |>
  opt_table_font(font = list("roboto"))
  # tab_style(style = cell_text(weight = "bold"),
  #           locations = cells_stub_grand_summary()) |>
  # tab_style(style = cell_text(weight = "bold"),
  #           locations = cells_column_labels()) |>
  # tab_style(style = cell_text(weight = "bold"),
  #           locations = cells_stubhead())

# Create a function to convert probabilities to odds:
probs.to.odds <- function (p) {
  stopifnot(p >= 0 & p <= 1)
  p/(1-p)}

# Check what the odds are for some probabilities:
probs.to.odds(c(0, 0.05, 0.5, 0.95, 0.99, 1))

# Make a sequence of probabilities from 0.001 to 0.999, incremented in steps of 0.001:
probabilities <- seq(0.001, 0.999, 0.001)

# Convert the probabilities to odds, using our function:
odds <- probs.to.odds(probabilities)

plot(x = probabilities, y = odds, ylim = c(0, 1)) # limit to odds under 1
plot(x = probabilities, y = odds, ylim = c(0, 100)) # limit to odds under 100
plot(x = probabilities, y = log(odds)) 
# instead of log(odds) could have used the built-in function logit: 
# plot(y = probabilities, x = logit(probabilities)) 
# check that the logit values are the same as the log(odds):
all.equal(logit(probabilities), log(odds))
#the logistic function:
inverselogit <- function (x) {1/(1+exp(-x))} 

# check that the logistic function really gets us back to probability:
all.equal(inverselogit(log(odds)), probabilities) 
titanic <- PASWR2::TITANIC3
titanic$sex <- factor(titanic$sex, levels=c("male","female"))

lm_titanic <- lm(survived ~ sex+age+pclass+sex:age, data = titanic)
car::Anova(lm_titanic, type = 3)

lm_titanic2 <- lm(survived ~ sex+age+pclass+sex:age+age:pclass+sex:pclass, data = titanic)
car::Anova(lm_titanic2, type = 3)

lm_titanic_plot <- marginaleffects::plot_predictions(lm_titanic2, by = c("age", "pclass", "sex"), 
                                  vcov = "HC3", points = 0.5, rug = FALSE) +
  geom_dots(alpha = 0.5, scale = .4,
    data = titanic, aes(x = age,y = survived, 
                        side = ifelse(survived == 1, "bottom", "top"))) +
  ylab("probability of survival") +
  scale_y_continuous(limits = c(-0.4, 1.2), breaks = c(0, 0.5, 1)) +
  scale_x_continuous(limits = c(0, 80)) +  
  scale_colour_manual(values = c( "#e6550d", "#4daf4a","#377eb8"),
                      breaks = c("1st", "2nd", "3rd"),
                      labels = c("1st", "2nd", "3rd")) +
  guides(colour = guide_legend(title = "class"),
         fill = guide_legend(title = "class")) 

glm_titanic1 <- glm(survived ~ sex*age*pclass,
                    family = binomial(link='logit'), data = titanic)

glm_titanic2 <- glm(survived ~ sex+age+pclass+sex:age+age:pclass,
                    family = binomial(link='logit'), data = titanic)

glm_titanic3 <- glm(survived ~ sex+age+pclass+sex:age, 
                    family = binomial(link='logit'), data = titanic)
car::Anova(glm_titanic3, type = 3)

glm_titanic4 <- glm(survived ~ sex+age+pclass, 
                    family = binomial(link='logit'), data = titanic)

glm_titanic5 <- glm(survived ~ sex+age+pclass+sex:age+age:pclass+sex:pclass, 
                    family = binomial(link='logit'), data = titanic)

AICcmodavg::aictab(list(glm_titanic1, glm_titanic2, glm_titanic3, glm_titanic4, glm_titanic5))

glm_titanic_plot <- marginaleffects::plot_predictions(glm_titanic5, by = c("age", "pclass", "sex"), 
                                  vcov = "HC3", points = 0.5, rug = FALSE) +
  geom_dots(alpha = 0.5, scale = .4,
    data = titanic, aes(x = age,y = survived, 
                        side = ifelse(survived == 1, "bottom", "top"))) +
  ylab("probability of survival") +
  scale_y_continuous(limits = c(-0.4, 1.2), breaks = c(0, 0.5, 1)) +
  scale_x_continuous(limits = c(0, 80)) +  
  scale_colour_manual(values = c( "#e6550d", "#4daf4a","#377eb8"),
                      breaks = c("1st", "2nd", "3rd"),
                      labels = c("1st", "2nd", "3rd")) + 
  guides(colour = guide_legend(title = "class"),
         fill = guide_legend(title = "class")) 

lm(survived ~ sex+age+pclass+sex:age+age:pclass+sex:pclass, data = titanic)
lm_titanic_plot 
glm(survived ~ sex+age+pclass+sex:age+age:pclass+sex:pclass, family = binomial(link = "logit"), data = titanic)
glm_titanic_plot 
#The placeholder there so that refined gets a box

ggplot() +
  geom_bar(
    data = pivoted,
    aes(
      x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, land_sea),
      fill = outcome,
      weight = antall
    ),
    position = "fill",
    linewidth = 0.3,
    width = 0.7,
    linejoin = 'mitre', lineend = 'square'
  ) +
  scale_fill_brewer(palette = "Blues", guide = guide_legend(title = "Outcome", order = 1)) +
  scale_color_brewer(palette = "Set1", direction = -1, breaks = c("r6", "0"), 
                     guide = guide_legend(title = paste0("EX rank<br>",
                                                        "<span style='font-size: 8pt'>",
                                                        "95% C.I.</span>"), 
                                          order = 2)) +
  geom_linerange(data = CIgs, #geom_pointrange
                  aes(
                    y = mean,
                    x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, land_sea),
                    colour = Ex_rank,
                    ymin = lower,
                    ymax = upper),
                 linewidth = 0.7) +
  scale_y_continuous(breaks = seq(0, 1, 0.2),
                     limits = c(0, 1),
                     expand = expansion(mult = c(0, 0.0))) +
  scale_linetype_discrete(name = NULL) +
  xlab("Treatment") +
  ylab("Proportion") +
  scale_x_discrete(guide = guide_axis_nested(title = NULL, type = "box", drop = FALSE)) +
  theme_guide(box = element_rect(colour = "grey92")) +
  annotate(geom = "richtext",
           size = 3.33,
           label = paste("**Factors**",
                         "<br>refined",
                         "<br>casting rank",
                         "<br>EX rank",
                         "<br>where"),
           fill = NA,
           label.color = NA,
           x = 12.4,
           y = 0,
           hjust = 0,
           vjust = 0.96) +
  coord_cartesian(expand = FALSE,
                  clip = 'off',
                  ylim = c(0, 1)) +
  geom_text(data = pivoted, aes(
      label =  paste0("n = ", antall),
      fill = outcome, #if dont have fill here, the labels will swap places! 
      y = antall / total,
      x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, land_sea)),
      position = position_fill(vjust = 0.5, reverse = FALSE),
      angle = 90,
      colour = "grey30",
      size = 3) +
 # annotation_custom(img, ymin = 0.8, ymax = 1, xmin = 12.8, xmax = 14.3) +
 # theme(axis.line.x.bottom = element_line(colour = "grey30")) 
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 9,
                                   vjust = 0.5,
                                   margin = margin(1.5, 0, 1.5, 0, "pt")),
        plot.margin = margin(c(10, 7, 7, 7), unit = "pt")) +
  theme(legend.title = element_markdown())

#for EX/no EX on land
ggplot() +
  geom_bar(
    data = EXpivoted,
    aes(
      x = Ex_rank,
      fill = outcome,
      #colour = Ex_rank,
      weight = antall
    ),
    position = "fill",
    linewidth = 0.3,
    width = 0.7,
    linejoin = 'mitre', lineend = 'square'
  ) +
  scale_fill_brewer(palette = "Blues", guide = guide_legend(title = "Outcome", order = 1)) +
  scale_color_brewer(palette = "Set1", direction = -1, breaks = c("r6", "0"), 
                     guide = guide_legend(title = paste0("EX rank<br>",
                                                        "<span style='font-size: 8pt'>",
                                                        "95% C.I.</span>"),
                                          order = 2)) +
  geom_linerange(data = EXCIgs,
                  aes(
                    y = mean,
                    x = Ex_rank,
                    ymin = lower,
                    ymax = upper,
                    #linetype = "95% C.I.",
                    colour = Ex_rank,
                  ),
                 linewidth = 0.8) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.2),
    limits = c(0, 1),
    expand = expansion(mult = c(0, 0.05))
  ) +
  scale_linetype_discrete(name = NULL) +
  xlab("EX rank") +
  ylab("Proportion") +
  coord_cartesian(expand = FALSE,
                  clip = 'off',
                  ylim = c(0, 1)) +
  geom_text(
    data = EXpivoted,
  aes(
      label =  paste0("n = ", antall),
      fill = outcome,
      y = 0.35, #y = antall / total,
      x = Ex_rank
    ),
    position = position_fill(vjust = 0.5, reverse = FALSE),
    angle = 0,
    colour = "grey30",
    size = 3,
  ) +
  geom_text(
    data = EXpivoted,
  aes(
      label =  paste0(round(antall/total*100, 1), "%"),
      fill = outcome,
      y = 0.35, 
      x = Ex_rank
    ),
    position = position_fill(vjust = 0.6, reverse = FALSE),
    angle = 0,
    colour = "grey30",
    size = 4,
  ) +
  #annotation_custom(img, ymin = 0.8, ymax = 1, xmin = 2.5, xmax = 3.1) +
  annotate(
    geom = "richtext",
    size = 3.2,
    label = paste("X<sup>2</sup> = ", #X<sup>2</sup><sub>1</sub> = 
                  round(EX_X2$statistic, 2), 
                  "<br>p = ", round(EX_X2$p.value, 3)),
    fill = NA,
    label.color = NA,
    x = 2.37,
    y = 0,
    hjust = 0,
    vjust = 1.15,
    colour = "grey30"
  ) +
  theme(legend.title = element_markdown(),
        plot.margin = margin(c(10, 7, 7, 7), unit = "pt"))
  
SeaEXcombinedProp <- combinedProp |>
  filter(land_sea == "sea") |>
  group_by(Ex_rank) |>
  summarise(normal = sum(normal),
            gs = sum(gs),
            total = sum(total))
SeaEXbinomCIgs <- binom::binom.wilson(SeaEXcombinedProp$gs, SeaEXcombinedProp$total)
SeaEXCIgs <- cbind(SeaEXcombinedProp, SeaEXbinomCIgs)
SeaEX_X2 <- prop.test(SeaEXcombinedProp$gs, SeaEXcombinedProp$total, alternative = "less")

SeaEXpivoted <- SeaEXCIgs |>
  dplyr::select(Ex_rank:total) |>
  pivot_longer(cols = c(normal, gs), names_to = "outcome", values_to = "antall")
SeaEXpivoted$outcome <- factor(SeaEXpivoted$outcome, levels = c("normal", "gs"))

ggplot() +
  geom_bar(
    data = SeaEXpivoted,
    aes(
      x = Ex_rank,
      fill = outcome,
      #colour = Ex_rank,
      weight = antall
    ),
    position = "fill",
    linewidth = 0.3,
    width = 0.7,
    linejoin = 'mitre', lineend = 'square'
  ) +
  scale_fill_brewer(palette = "Blues", guide = guide_legend(title = "Outcome", order = 1)) +
  scale_color_brewer(palette = "Set1", direction = -1, breaks = c("r6", "0"),
                     guide = guide_legend(title = paste0("EX rank<br>",
                                                        "<span style='font-size: 8pt'>",
                                                        "95% C.I.</span>"),
                                          order = 2)) + 
  geom_linerange(data = SeaEXCIgs,
                  aes(
                    y = mean,
                    x = Ex_rank,
                    ymin = lower,
                    ymax = upper,
                    #linetype = "95% C.I.",
                    colour = Ex_rank,
                  ),
                 linewidth = 0.8) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.2),
    limits = c(0, 1),
    expand = expansion(mult = c(0, 0.05))
  ) +
  scale_linetype_discrete(name = NULL) +
  xlab("EX rank") +
  ylab("Proportion") +
  coord_cartesian(expand = FALSE,
                  clip = 'off',
                  ylim = c(0, 1)) +
  geom_text(
    data = SeaEXpivoted,
  aes(
      label =  paste0("n = ", antall),
      fill = outcome,
      y = 0.35, #y = antall / total,
      x = Ex_rank
    ),
    position = position_fill(vjust = 0.5, reverse = FALSE),
    angle = 0,
    colour = "grey30",
    size = 3,
  ) +
  geom_text(
    data = SeaEXpivoted,
  aes(
      label =  paste0(round(antall/total*100, 1), "%"),
      fill = outcome,
      y = 0.35, 
      x = Ex_rank
    ),
    position = position_fill(vjust = 0.6, reverse = FALSE),
    angle = 0,
    colour = "grey30",
    size = 4,
  ) +
  #annotation_custom(img, ymin = 0.8, ymax = 1, xmin = 2.5, xmax = 3.1) +
  annotate(
    geom = "richtext",
    size = 3.2,
    label = paste("X<sup>2</sup> = ", #X<sup>2</sup><sub>1</sub> = 
                  round(SeaEX_X2$statistic, 2), 
                  "<br>p = ", round(SeaEX_X2$p.value, 3)),
    fill = NA,
    label.color = NA,
    x = 2.37,
    y = 0,
    hjust = 0,
    vjust = 1.15,
    colour = "grey30"
  ) +
  theme(legend.title = element_markdown(),
        plot.margin = margin(c(10, 7, 7, 7), unit = "pt"))

#combine all groups, also at sea. But might be confounded because one more group at sea with EX than without so might be effect of workshop instead

##slå sammen de med ex og de uten, INKLUDERT land og lag nye CIs
AllEXcombinedProp <- combinedProp |>
  group_by(Ex_rank) |>
  summarise(normal = sum(normal),
            gs = sum(gs),
            total = sum(total))
AllEXbinomCIgs <- binom::binom.wilson(AllEXcombinedProp$gs, AllEXcombinedProp$total)
AllEXCIgs <- cbind(AllEXcombinedProp, AllEXbinomCIgs)
AllEX_X2 <- prop.test(AllEXcombinedProp$gs, AllEXcombinedProp$total, alternative = "less")

AllEXpivoted <- AllEXCIgs |>
  dplyr::select(Ex_rank:total) |>
  pivot_longer(cols = c(normal, gs), names_to = "outcome", values_to = "antall")
AllEXpivoted$outcome <- factor(AllEXpivoted$outcome, levels = c("normal", "gs"))

ggplot() +
  geom_bar(
    data = AllEXpivoted,
    aes(
      x = Ex_rank,
      fill = outcome,
      #colour = Ex_rank,
      weight = antall
    ),
    position = "fill",
    linewidth = 0.3,
    width = 0.7,
    linejoin = 'mitre', lineend = 'square'
  ) +
  scale_fill_brewer(palette = "Blues", guide = guide_legend(title = "Outcome", order = 1)) +
  scale_color_brewer(palette = "Set1", direction = -1, breaks = c("r6", "0"),
                     guide = guide_legend(title = paste0("EX rank<br>",
                                                        "<span style='font-size: 8pt'>",
                                                        "95% C.I.</span>"),
                                          order = 2)) +
  geom_linerange(data = AllEXCIgs,
                  aes(
                    y = mean,
                    x = Ex_rank,
                    ymin = lower,
                    ymax = upper,
                    #linetype = "95% C.I.",
                    colour = Ex_rank,
                  ),
                 linewidth = 0.8) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.2),
    limits = c(0, 1),
    expand = expansion(mult = c(0, 0.05))
  ) +
  scale_linetype_discrete(name = NULL) +
  xlab("EX rank") +
  ylab("Proportion") +
  coord_cartesian(expand = FALSE,
                  clip = 'off',
                  ylim = c(0, 1)) +
  geom_text(
    data = AllEXpivoted,
  aes(
      label =  paste0("n = ", antall),
      fill = outcome,
      y = 0.35, #y = antall / total,
      x = Ex_rank
    ),
    position = position_fill(vjust = 0.5, reverse = FALSE),
    angle = 0,
    colour = "grey30",
    size = 3,
  ) +
  geom_text(
    data = AllEXpivoted,
  aes(
      label =  paste0(round(antall/total*100, 1), "%"),
      fill = outcome,
      y = 0.35, 
      x = Ex_rank
    ),
    position = position_fill(vjust = 0.6, reverse = FALSE),
    angle = 0,
    colour = "grey30",
    size = 4,
  ) +
  annotate(
    geom = "richtext",
    size = 3.2,
    label = paste("X<sup>2</sup> = ", #X<sup>2</sup><sub>1</sub> = 
                  round(AllEX_X2$statistic, 2), 
                  "<br>p = ", round(AllEX_X2$p.value, 3)),
    fill = NA,
    label.color = NA,
    x = 2.37,
    y = 0,
    hjust = 0,
    vjust = 1.15,
    colour = "grey30"
  ) +
  theme(legend.title = element_markdown(),
        plot.margin = margin(c(10, 7, 7, 7), unit = "pt"))

# successrate------------------


#statistics

# create 'result' vector
# repeat 1s and 0s the number of times given in the respective 'count' column
# result <- rep(rep(c(1,0), nrow(additional)), unlist(additional[ , c("nGs", "nNormal")]))
# 
# # repeat each row in df the number of times given by the sum of 'count' columns
# Longadditional <- data.frame(additional[rep(1:nrow(additional), rowSums(additional[ , c("nGs", "nNormal")]) ), c("sailor_equip", "refined", "crafting_rank", "Ex_rank")], result)
# Longadditional$outcome <- factor(Longadditional$outcome, levels = c("normal", "gs"))

# invert_table <- function(x){
#       y <- x[rep(rownames(x),x$Freq),1:(ncol(x)-1)]
#       rownames(y) <- c(1:nrow(y))
#       return(y)}
# invert_table(additional)

#setNames(tidyr::uncount(as.data.frame(tab), Freq), c("Type", "Origin"))

normalAdditional <- additional |>
  uncount(nNormal) |> mutate(outcome = "normal") |> select(id:pp, comm:placeholder)

gsAdditional <- additional |>
  uncount(nGs) |> mutate(outcome = "gs") |> select(id:pp, comm:placeholder)

Longadditional <- bind_rows(gsAdditional, normalAdditional)

Longminions <- bind_rows(minions, Longadditional)
Longminions$land_sea <- factor(Longminions$land_sea, levels = c("land", "sea"))
Longminions$refined <- factor(Longminions$refined)
Longminions$crafting_rank <- factor(Longminions$crafting_rank, levels = c("r3", "r20"))
Longminions$EX <- factor(Longminions$Ex_rank)
Longminions$where <- Longminions$land_sea
Longminions$outcome <- factor(Longminions$outcome, levels = c("normal", "gs"))
attr(Longminions$land_sea, "label") <- "where "

options(contrasts=c("contr.sum","contr.poly"))
gt_glmTable <- function(df, event, title, keep = FALSE, ...) {
  # arguments: model, failure/great success, title, keep (p values for each level)
  tbl_regression(df,
                 exponentiate = TRUE,
                 add_estimate_to_reference_rows = FALSE) |>
    add_n(location = c('level')) |>
    add_nevent(location = c('level')) |>
    add_global_p(keep = keep, type = "III") |>
    # this will update the em-dash in the CI row to Ref.
    # modify_table_styling(
    #   columns = conf.low,
    #   rows = reference_row %in% TRUE,
    #   missing_symbol = "Reference"
    # ) |>
    # adding event rate
    modify_table_body( ~ .x |>
                         dplyr::mutate(
                           stat_nevent_rate =
                             ifelse(!is.na(stat_nevent), paste0(
                               style_sigfig(stat_nevent / stat_n, scale = 100, digits = 3), "%"
                             ), NA),
                           .after = stat_nevent)) |>
    # merge the colums into a single column
    modify_column_merge(pattern = "{stat_nevent} / {stat_n} ({stat_nevent_rate})", 
                        rows = !is.na(stat_nevent)) |>
    modify_table_styling(column = conf.low,
                         rows = !is.na(estimate),
                         cols_merge_pattern = "{conf.low}–{conf.high}") |>
    modify_header(stat_nevent = glue::glue("{event} Rate")) |>
    bold_labels() |>
    modify_header(label = "Variable", #to remove bold in headers
                  estimate = "OR", 
                  conf.low = "95% CI", 
                  p.value = "p-value") |> 
    modify_column_alignment(columns = stat_nevent, align = "right") |>
    as_gt() |>
    tab_header(title = title) |>
    opt_stylize(style = 3, color = "blue") |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
    opt_horizontal_padding(scale = 1.5) |>
    tab_options(data_row.padding = px(2),
                summary_row.padding = px(3),
                grand_summary_row.padding = px(3),
                row_group.padding = px(4),
                footnotes.padding = px(2),
                source_notes.padding = px(2)) |>
    tab_style(locations = cells_title (), style = cell_borders(style = 'hidden')) |>
    cols_align_decimal(columns = p.value,
                       dec_mark = ".",
                       locale = NULL) |>
    tab_style(style = cell_text(align = "center"),
              locations = cells_column_labels(columns = stat_nevent)) |>
    tab_style(style = cell_text(size = px(15)),
              locations = cells_body()) |>
    tab_style(style = cell_text(size = px(12)),
              locations = cells_footnotes()) |>
    tab_style(style = cell_text(size = px(12), align = "left"),
              locations = cells_source_notes()) |>
    opt_table_font(font = list("roboto"))
}

gt_anovaTable <- function(df, title){
gt(df, rowname_col = "Effect") |>
  tab_stubhead(label = "Effect") |>
  tab_footnote(footnote = "partial eta-squared (effect size).", 
               locations = cells_column_labels(columns = pes)) |>
  tab_style(style = cell_text(size = px(13)),
            locations = cells_footnotes()) |>
  tab_header(title = glue::glue("{title}")) |>
  tab_style(style = cell_text(color = "#777"), locations = cells_title ()) |>
  cols_align(align = 'right', columns = 3:5) |> 
  cols_align(align = 'center', columns = 2) |>
  opt_stylize(style = 3, color = "blue") |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>  
  tab_options(data_row.padding = px(2),
              summary_row.padding = px(3),
              grand_summary_row.padding = px(3),
              row_group.padding = px(4),
              footnotes.padding = px(2))|>
  opt_horizontal_padding(scale = 1.5) |>
  tab_style(locations = cells_title (), style = cell_borders(style = 'hidden')) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_body()) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_stub())  |>
  # tab_style(locations =  cells_column_labels(),
  #           style = list(cell_text(weight = 'bold'))) |>
  # tab_style(locations =  cells_stubhead(),
  #           style = list(cell_text(weight = 'bold'))) |>
  text_replace(pattern = "_",
               replacement = " ",
               locations = cells_stub())|>
  text_replace(pattern = ":",
               replacement = " x ",
               locations = cells_stub()) |>
  text_replace(pattern = "\\.", 
               replacement = "0.",
               locations = cells_body(columns = c(4,5))) |>
  cols_label(p.value = "p") |>
  opt_table_font(font = list("roboto"))
}
options(contrasts=c("contr.sum","contr.poly"))

glm1 <- glm(outcome ~ sailor_equip+Ex_rank,
            family = binomial(link='logit'), data = Longminions)

summary(glm1)
car::Anova(glm1, type = 3)
anova(glm1)

# glm1 |> broom.helpers::tidy_plus_plus(
#     exponentiate = TRUE,
#     add_reference_rows = FALSE,
#     categorical_terms_pattern = "{level} / {reference_level}",
#     add_n = TRUE
#   ) |> gt()

gt_glmTable(glm1, "Great Success", glm1$formula)

options(contrasts=c("contr.treatment","contr.poly")) #dummy coding

glmNULL <- glm(outcome~1,  family = binomial(link='logit'), data = Longminions)

glmFULL <- glm(outcome ~ EX*where*crafting_rank*refined,
            family = binomial(link='logit'), data = Longminions)

glmADDITIVE <- glm(outcome ~ EX+where+crafting_rank+refined,
            family = binomial(link='logit'), data = Longminions)
car::Anova(glmADDITIVE, type = 3)
gt_glmTable(glmADDITIVE, "Great Success", glm1$formula)

glm_no_refined_full <- glm(outcome ~ EX*where*crafting_rank,
            family = binomial(link='logit'), data = Longminions)

glm_no_refined_additive <- glm(outcome ~ EX+where+crafting_rank,
            family = binomial(link='logit'), data = Longminions)

glm_no_refined_no_casting_full <- glm(outcome ~ EX*where,
            family = binomial(link='logit'), data = Longminions)


glm_no_refined_no_casting_additive <- glm(outcome ~ EX+where,
            family = binomial(link='logit'), data = Longminions)
car::Anova(glm_no_refined_no_casting_additive, type = 3)
summary(glm_no_refined_no_casting_additive)
#plot(glm_no_refined_no_casting_additive)

#kanskje endre navn og bruke land_sea istedet, ettersom det er det som er effekten

predictionplotGlm <- 
marginaleffects::plot_predictions(glm_no_refined_no_casting_additive,
                                  condition = c("where", "EX"),
                                  #coef_rename = TRUE
                                  ) +
  ylab("Predicted probability of great success") +
  xlab("On land (with Meister title and apprentice) or \nat sea (with workshop and sailor equip.)") +
  scale_colour_manual(values = c("#e41a1c", "#377eb8"), 
                      labels = c("r6","0"),
                      breaks = c("r6", "0")) +
  scale_y_continuous(breaks = seq(0.44, 0.56, 0.02), 
                     limits = c(0.44, 0.56)) +
    guides(colour = guide_legend(title = "EX rank", order = 1),
           alpha = guide_legend(Title = NULL, 
                                override.aes = list(alpha = 1),
                                order = 2)) 

glm_onlySE <- glm(outcome ~ where,
            family = binomial(link='logit'), data = Longminions)
car::Anova(glm_onlySE, type = 3)

glm_onlyEX <- glm(outcome ~ EX,
            family = binomial(link='logit'), data = Longminions)
car::Anova(glm_onlyEX, type = 3)

AIC(glmFULL, glmADDITIVE, glm_no_refined_full, glm_no_refined_additive, glm_no_refined_no_casting_full, glm_no_refined_no_casting_additive, glm_onlyEX, glm_onlySE, glmNULL)

models_glm <-list(glmFULL, glm_no_refined_full, glmADDITIVE, glm_no_refined_no_casting_full, glm_no_refined_additive, glm_no_refined_no_casting_additive, glm_onlySE, glm_onlyEX, glmNULL)

modnames_glm <- paste0("(", 1:length(models_glm), ")")

modelSummaryGLM <- modelsummary(models_glm,   
             estimate  = "{estimate}\n[{conf.low}, {conf.high}]\np = {p.value}",  
             statistic = NULL,
             fmt = fmt_decimal(digits = 2, pdigits = 3),
             #add_rows = emptyrow,
             #coef_omit = "Intercept",
             gof_omit = "Log.Lik.|RMSE|F",
             exponentiate = TRUE,
             output = "gt",
             include_reference = TRUE,
             coef_rename = TRUE,
             #coef_rename = coef_rename
             # coef_rename = c("sailor_Equip" = "Sailor Equip. ", 
             #                 "refined" = "Refined ",
             #                 "crafting_rank" = "Casting Rank ",
             #                 "Ex_rank" = "EX Rank "),
             #rowname_col = c(1) #for some reason the GT options are not honored
             #rownames_to_stub = TRUE
             ) |>   
tab_stubhead(label = "Parameters") |>
    cols_label(c(1) ~ md("**Parameters**")) |>
      # tab_spanner(label = md("**Chosen**"), 
      #         columns = "(6)",
      #         id = "d") |>
  tab_spanner(label = md("**Candidate Models (logistic regression)**"), 
              columns = c(2:10),
              id = "c") |>
   # tab_spanner(label = md("**Null Model**"), columns = c(2)) |>
    tab_footnote(footnote = "× indicates interaction", 
               locations = cells_column_labels(columns=c(1))) |>
    tab_footnote(footnote = "Values shown are the Odds Ratios and [95% confidence intervals]. Contrasts: dummy coding.", 
               locations = cells_column_spanners(spanners = "c")) |>
  opt_stylize(style = 3, color = "blue") |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
  tab_options(
      data_row.padding = px(2),
      summary_row.padding = px(3),
      grand_summary_row.padding = px(3),
      row_group.padding = px(4),
      footnotes.padding = px(2),
      quarto.disable_processing = TRUE
    ) |>
  cols_width(starts_with("(") ~ pct(9)) |>
  tab_style(
    style = cell_text(size = px(13)),
    locations = cells_body()
    ) |>
   tab_style(
    style = cell_text(size = px(11)),
    locations = cells_footnotes()
    ) |>
   tab_style(
    locations =  cells_body(columns = c(1), rows = c(1:16)),
    style = list(cell_fill(color = 'grey85'),
                 cell_text(weight = 'bold'))
  ) |>
   sub_values(
  values = "-",
  replacement = "—",
  escape = TRUE
)  |>
  tab_style(
    locations =  cells_body(columns = c(2:10)),
    style = cell_text(whitespace = 'pre')
  )|>
   tab_style(
    cell_fill(color = "#deebf7"),
    location = cells_body(columns = "(6)")
  ) |>
    opt_table_font(font = list("roboto"))

modelSummaryGLM |> tab_info()

AICcmodavg::aictab(list(glmFULL, glm_no_refined_full, glmADDITIVE, glm_no_refined_no_casting_full, glm_no_refined_additive, glm_no_refined_no_casting_additive, glm_onlySE, glm_onlyEX, glmNULL))

AICglm <- AICcmodavg::aictab(models_glm, 
                            modnames = modnames_glm) |>
  gt() |>
    opt_stylize(style = 3, color = "blue") |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
    tab_options(
      data_row.padding = px(2),
      summary_row.padding = px(3),
      grand_summary_row.padding = px(3),
      row_group.padding = px(4),
      footnotes.padding = px(2),
      quarto.disable_processing = TRUE
    ) |> 
  opt_horizontal_padding(scale = 1.5) |>
    cols_align(align = 'center', columns = 1) |> 
  cols_label(
    Modnames = "Model",
    Delta_AICc = "ΔAICc"
  ) |>
  fmt_number(decimals = 2, columns = c(3, 4, 5, 6, 8), use_seps = FALSE) |>
   cols_hide(columns = c(5:8)) |>
   tab_style(
    cell_fill(color = "#deebf7"),
    location = cells_body(rows =  c(2))
  ) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_body()) |>
  opt_table_font(font = list("roboto"))

glmTable_succ <- gt_glmTable(glm_no_refined_no_casting_additive, "Great Success", NULL) |>
      tab_footnote(footnote = "Contrasts: dummy-coding", 
               locations = cells_column_labels(columns = "estimate")) |>
     tab_style(
       style = cell_text(size = px(12)),
       locations = cells_footnotes()) |>
     tab_options(quarto.disable_processing = TRUE)

car::Anova(glm_no_refined_no_casting_additive, type = 3)
anova(glm_no_refined_no_casting_additive)

gt_glmTable(glm1, "Great Success", NULL) 
car::Anova(glm1) 

boot.pval::boot_summary(glm_no_refined_no_casting_additive, type = "perc", method = "case", coef ="exp")
arm::binnedplot(predict(glm_no_refined_no_casting_additive, type = "response"),
                residuals(glm_no_refined_no_casting_additive, type = "response"))

#PP-----------------------

options(contrasts=c("contr.sum","contr.poly"))

lm1pp <- lm(pp ~ sailor_equip*refined*crafting_rank+Ex_rank, 
                  data = subset(minions, durability >= 100))
car::Anova(lm1pp, type = 3)
summary(lm1pp)

lm2pp <- lm(pp ~ sailor_equip+refined+crafting_rank+Ex_rank, 
                  data = subset(minions, durability >= 100))
car::Anova(lm2pp, type = 3)

lm3pp <- lm(pp ~ sailor_equip+crafting_rank+Ex_rank, 
                  data = subset(minions, durability >= 100))

lm4pp <- lm(pp ~ sailor_equip+Ex_rank, 
                  data = subset(minions, durability >= 100))

lm5pp <- lm(pp ~ sailor_equip, 
                  data = subset(minions, durability >= 100))

lm6pp <- lm(pp ~ 1, 
                  data = subset(minions, durability >= 100))

models_pp <-list(lm1pp, lm2pp, lm3pp, lm4pp, lm5pp, lm6pp)

modnames_pp <- paste0("(", 1:length(models_pp), ")")

AICpp <- AICcmodavg::aictab(models_pp, 
                            modnames = modnames_pp) |>
  gt() |>
    opt_stylize(style = 3, color = "blue") |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
    tab_options(
      data_row.padding = px(2),
      summary_row.padding = px(3),
      grand_summary_row.padding = px(3),
      row_group.padding = px(4),
      footnotes.padding = px(2),
      quarto.disable_processing = TRUE) |>
  opt_horizontal_padding(scale = 1.5) |>
    cols_align(align = 'center', columns = 1) |> 
  cols_label(
    Modnames = "Model",
    Delta_AICc = "ΔAICc") |>
  fmt_number(decimals = 2, columns = c(3, 4,5,6,8), use_seps = FALSE) |>
   cols_hide(columns = c(5:8)) |>
   tab_style(
    cell_fill(color = "#deebf7"),
    location = cells_body(rows =  c(1))) |>
    tab_style(style = cell_text(size = px(15)),
            locations = cells_body()) |>
    opt_table_font(font = list("roboto"))

anova(lm5pp, lm6pp)

anovaTable1pp <- aov_car(pp ~ sailor_equip + Error(id),
                  data = subset(minions, pp >= 100))
anovaTable2pp <- nice(anovaTable1pp, MSE = FALSE, sig_symbols = rep("", 4), es = "pes")

anovaTable_pp <- gt_anovaTable(anovaTable2pp, NULL)
  # text_replace(
  #   pattern = "\\.",
  #   replacement = "0.",
  #   locations = cells_body(columns = c(4,5)))

t.test(pp ~ sailor_equip, data=minions, var.equal = FALSE)

predictionplotlmpp <- 
marginaleffects::plot_predictions(lm5pp,
                                  by = c("sailor_equip"),
                                  vcov = "HC3"
                                  #coef_rename = TRUE
                                  ) +
  ylab("Predicted piercing power") +
  xlab("Sailor Equipment") + geom_pointrange(aes(y= estimate, x = sailor_equip, 
                                                 ymin = conf.low, ymax = conf.high,
                                                 colour = sailor_equip),
                                        show.legend = FALSE) +
    scale_colour_manual(values = c("#377eb8", "#e41a1c"))

modelSummaryPp <- modelsummary(models_pp,   
             estimate  = "{estimate}\n[{conf.low}, {conf.high}]\np = {p.value}",
             statistic = NULL,
             fmt = fmt_decimal(digits = 2, pdigits = 3),
             #add_rows = emptyrow,
             #coef_omit = "Intercept",
             gof_omit = "Log.Lik.|RMSE|F",
             exponentiate = FALSE,
             output = "gt",
             include_reference = TRUE,
             vcov = "robust",
             #shape =  ~ term  ~ model,
             #coef_rename = TRUE,
             coef_rename = coef_rename,
             # coef_rename = c("sailor_Equip" = "Sailor Equip. ", 
             #                 "refined" = "Refined ",
             #                 "crafting_rank" = "Casting Rank ",
             #                 "Ex_rank" = "EX Rank "),
             #rowname_col = c(1) #for some reason the GT options are not honored
             #rownames_to_stub = TRUE
             ) |>   
tab_stubhead(label = "Parameters") |>
  cols_label(c(1) ~ md("**Parameters**")) |>
  # tab_spanner(label = md("**Chosen**"), 
  #             columns = "(5)",
  #             id = "d") |>
  tab_spanner(label = md("**Candidate Models**"), 
              columns = c(2:7),
              id = "c") |>
  tab_footnote(footnote = "× indicates interaction", 
               locations = cells_column_labels(columns=c(1))) |>
  tab_footnote(footnote = "Contrasts: effects coding. Variance-covariance matrix: HC3 (standard errors robust to heteroskedasticity)", 
               locations = cells_column_spanners(spanners = "c")) |>
  opt_stylize(style = 3, color = "blue") |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
  tab_options(
      data_row.padding = px(2),
      summary_row.padding = px(3),
      grand_summary_row.padding = px(3),
      row_group.padding = px(4),
      footnotes.padding = px(2),
      quarto.disable_processing = TRUE
    ) |>
  cols_width(starts_with("(") ~ pct(12)) |>
  tab_style(
    style = cell_text(size = px(13)),
    locations = cells_body()
    ) |>
   tab_style(
    style = cell_text(size = px(11)),
    locations = cells_footnotes()
    ) |>
   tab_style(
    locations =  cells_body(columns = c(1), rows = c(1:9)),
    style = list(cell_fill(color = 'grey85'),
                 cell_text(weight = 'bold'))
  ) |>
    tab_style(
    locations =  cells_body(columns = c(2:7)),
    style = cell_text(whitespace = 'pre')
  ) |>
   sub_values(
  values = "-",
  replacement = "—",
  escape = TRUE
  ) |>
  text_replace(
    pattern = "_",
    replacement = " ",
    locations = cells_body()
  ) |>
  text_replace(
    pattern = ":",
    replacement = " × ",
    locations = cells_body()
  ) |>
   tab_style(
    cell_fill(color = "#deebf7"),
    location = cells_body(columns = "(5)")
  ) |>
    opt_table_font(font = list("roboto"))
  
  #   text_replace(
  #   pattern = "/^$/",
  #   replacement = "---",
  #   locations = cells_body()
  # ) |>
  #fmt_markdown() #NB! fmt_markdown fucker opp i quarto! så istedet for <br> bruk \n og whitespace-pre CSS
#Dura----------------------

minions$refined2 <- minions$refined
levels(minions$refined2) <- c("unrefined", "refined")

options(contrasts=c("contr.sum","contr.poly"))

lm1 <- lm(durability ~ sailor_equip*refined*crafting_rank+Ex_rank, 
                  data = subset(minions, durability >= 100),
                  )
summary(lm1)
ggResidpanel::resid_panel(lm1, plots = c("resid", "qq", "hist", "yvp"), nrow = 2)
marginaleffects::plot_predictions(lm1,
                                  condition = c("sailor_equip", "refined", "crafting_rank"),
                                  #coef_rename = TRUE
                                  )

#plot(lm1)
jtools::effect_plot(model = lm1, pred = crafting_rank, interval = TRUE, partial.residuals = TRUE, jitter = .2)

lmtest::bptest(lm1)

lmtest::coeftest(lm1, vcov = vcovHC(lm1, type = "HC3"))

lm2 <- lm(durability ~ sailor_equip+refined+crafting_rank+Ex_rank,
                  data = subset(minions, durability >= 100)
                  )
ggResidpanel::resid_panel(lm2, plots = c("resid", "qq", "hist", "yvp"), nrow = 1)
jtools::effect_plot(model = lm1, pred = refined, interval = TRUE, partial.residuals = TRUE, jitter = .2)

lm3 <- lm(durability ~ sailor_equip*refined*crafting_rank,
                  data = subset(minions, durability >= 100)
                  )
#to get level names unrefined and refined instead of changing many places.
lm3b <- lm(durability ~ sailor_equip*refined2*crafting_rank,
                  data = subset(minions, durability >= 100)
                  )

lm4 <- lm(durability ~ sailor_equip+refined+crafting_rank,
                  data = subset(minions, durability >= 100)
                  )

lm5 <- lm(durability ~ sailor_equip+refined + crafting_rank + sailor_equip:refined + sailor_equip:crafting_rank + refined:crafting_rank,
                  data = subset(minions, durability >= 100)
                  )

lm6 <- lm(durability ~ 1,
                  data = subset(minions, durability >= 100)
                  )
models_dura <-list(lm1, lm2, lm3, lm4, lm5, lm6)

modelSummaryDura <- modelsummary(models_dura,   
             estimate  = "{estimate}\n[{conf.low}, {conf.high}]\np = {p.value}",  
             statistic = NULL,
             fmt = fmt_decimal(digits = 2, pdigits = 3),
             #add_rows = emptyrow,
             #coef_omit = "Intercept",
             gof_omit = "Log.Lik.|RMSE|F",
             exponentiate = FALSE,
             output = "gt",
             include_reference = TRUE,
             vcov = "robust",
             shape =  ~ term  ~ model,
             #coef_rename = TRUE
             coef_rename = coef_rename,
             # coef_rename = c("sailor_Equip" = "Sailor Equip. ", 
             #                 "refined" = "Refined ",
             #                 "crafting_rank" = "Casting Rank ",
             #                 "Ex_rank" = "EX Rank "),
             #rowname_col = c(1) #for some reason the GT options are not honored
             #rownames_to_stub = TRUE
             ) |>   
tab_stubhead(label = "Parameters") |>
    cols_label(c(1) ~ md("**Parameters**")) |>
    # tab_spanner(label = md("**Chosen**"), 
    #           columns = "(3)",
    #           id = "d") |>
  tab_spanner(label = md("**Candidate Models**"), 
              columns = c(2:7),
              id = "c") |>
   # tab_spanner(label = md("**Null Model**"), columns = c(2)) |>
  tab_footnote(footnote = "× indicates interaction", 
               locations = cells_column_labels(columns=c(1))) |>
  tab_footnote(footnote = "Contrasts: effects coding. Variance-covariance matrix: HC3 (standard errors robust to heteroskedasticity)", 
               locations = cells_column_spanners(spanners = "c")) |>
  opt_stylize(style = 3, color = "blue") |> 
  tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
  tab_options(
      data_row.padding = px(2),
      summary_row.padding = px(3),
      grand_summary_row.padding = px(3),
      row_group.padding = px(4),
      footnotes.padding = px(2),
      quarto.disable_processing = TRUE
    ) |>
  #cols_width(starts_with("(") ~ pct(9)) |>
  tab_style(
    style = cell_text(size = px(13)),
    locations = cells_body()
    ) |>
   tab_style(
    style = cell_text(size = px(11)),
    locations = cells_footnotes()
    ) |>
   tab_style(
    locations =  cells_body(columns = c(1), rows = c(1:9)),
    style = list(cell_fill(color = 'grey85'),
                 cell_text(weight = 'bold'))
  ) |>
    tab_style(
    locations =  cells_body(columns = c(2:7)),
    style = cell_text(whitespace = 'pre')
  ) |>
   sub_values(
  values = "-",
  replacement = "—",
  escape = TRUE
  ) |>
  text_replace(
    pattern = "_",
    replacement = " ",
    locations = cells_body()
  ) |>
  text_replace(
    pattern = ":",
    replacement = " × ",
    locations = cells_body()
  ) |>
   tab_style(
    cell_fill(color = "#deebf7"),
    location = cells_body(columns = "(3)")
  )  |>
    opt_table_font(font = list("roboto"))


AICcmodavg::aictab(list(lm1, lm2, lm3, lm4, lm5, lm6))

modnames_dura <- paste0("(", 1:length(models_dura), ")")

AICdura <- AICcmodavg::aictab(models_dura, 
                            modnames = modnames_dura) |>
  gt() |>
    opt_stylize(style = 3, color = "blue") |>  
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_body()) |>
    tab_style(style = cell_borders(sides = c("top"),  color = NULL, weight = px(0.5)),
            locations = cells_stub()) |>
    tab_options(
      data_row.padding = px(2),
      summary_row.padding = px(3),
      grand_summary_row.padding = px(3),
      row_group.padding = px(4),
      footnotes.padding = px(2),
      quarto.disable_processing = TRUE
    ) |> 
  opt_horizontal_padding(scale = 1.5) |>
    cols_align(align = 'center', columns = 1) |> 
  cols_label(
    Modnames = "Model",
    Delta_AICc = "ΔAICc"
  ) |>
  fmt_number(decimals = 2, columns = c(3, 4,5,6,8), use_seps = FALSE) |>
  cols_hide(columns = c(5:8)) |>
  tab_style(
    cell_fill(color = "#deebf7"),
    location = cells_body(rows =  c(1))
  ) |>
  tab_style(style = cell_text(size = px(15)),
            locations = cells_body()) |>
  opt_table_font(font = list("roboto"))

predictionplotlm <- 
marginaleffects::plot_predictions(lm3b,
                                  condition = c("crafting_rank", "sailor_equip", "refined2"),
                                  vcov = "HC3"
                                  #coef_rename = TRUE
                                  ) +
  ylab("Predicted duralibity") +
  xlab("Casting rank") +
  scale_colour_manual(values = c( "#e41a1c", "#377eb8"),
                      breaks = c("r5", "(land)"),
                      labels = c( "r5", "(on land)"))  +
  guides(colour = guide_legend(title = "Sailor Equip."))+
  scale_y_continuous(breaks = seq(115, 140, 5), 
                     limits = c(114, 140)) 
  # labs(subtitle = "Casting rank") + 
  # theme(plot.subtitle = element_textbox_simple(halign = 0.5,
  #                                              fill = "gray80",
  #                                              padding = margin(3,3,3,3))) 

anovaTable1 <- aov_car(durability ~ sailor_equip*refined*crafting_rank + Error(id),
                  data = subset(minions, durability >= 100))
anovaTable2 <- nice(anovaTable1, MSE = FALSE, sig_symbols = rep("", 4), es = "pes")

#Anova Table (Type III tests) of the effect of various factors on the durability

anovaTable_dura <- gt_anovaTable(anovaTable2, NULL) |>   
  # text_replace(
  #   pattern = "\\.",
  #   replacement = "0.",
  #   locations = cells_body(columns = c(4,5))
  # ) |>
  #tab_options(table.width = pct(50)) |>
  tab_footnote(footnote = "x indicates interaction.", 
               locations = cells_stubhead()) |>
  tab_options(footnotes.multiline = FALSE,
  footnotes.sep = " ",) |>
  cols_width(c(df) ~ pct(16))

glmTable_succ
predictionplotGlm
ppmean <- round(mean(minions$pp, na.rm = T), 2)

#GeomHline$default_aes

minions |>
  ggplot(aes(y = pp, x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, sailor_equip))) +
  # geom_boxplot(,#aes(fill = Ex_rank),
  #              outliers = F,
  #              staplewidth = 0.7,
  #              #linejoin = 'mitre', lineend = 'square'
  #              ) +
  geom_count(alpha = 0.6, shape = 16, colour = "grey70") +
  scale_size_area() +
  ylab("Piercing Power") +
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 9,
                                   vjust = 0.5, 
                                   margin = margin(1.5,0,1.5,0, "pt")),
        legend.key.spacing.y = unit(0.5, "pt")) +
  scale_x_discrete(guide = guide_axis_nested(
    title = NULL,
    type = "box",
    drop = FALSE,
    #levels_text =
  )) +
  theme_guide(box = element_rect(colour = "grey92")) +
  annotate(
    geom = "richtext",
    size = 3.33,
    label = paste("**Factors**", "<br>refined", "<br>casting rank", 
                  "<br>EX rank", "<br>sailor equipment"),
    fill = NA,
    label.color = NA,
    x = 12.4,
    y = 152.15,
    hjust = 0,
    vjust = 1
  ) + 
  scale_y_continuous(breaks = c(152, 156, 160, 164, 168),
                         expand = expansion(mult = c(0, 0.1))) +
  coord_cartesian(expand = FALSE,
                  clip = 'off',
                  ylim = c(152, 168),
                  xlim = c(0.7, 12.3)) +
  geom_hline(aes(yintercept = mean(pp, na.rm = T), 
                 alpha = paste("overall mean\n=", ppmean)),
              linetype = 2,
              linewidth = 0.25,
              colour = "grey30",
              #alpha = 1
              ) +
  geom_linerange(data = pldatPp, aes(
    y = mean,
    x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, sailor_equip),
    ymin = low,
    ymax = upp,
    linetype = "mean and \n95% C.I.\nbootstrap",
    colour = sailor_equip)
    ) +
    scale_linetype_discrete(name = NULL) + #name = NULL
  stat_summary(
    fun = mean,
    aes(shape = "mean and \n95% C.I.\nbootstrap", fill = sailor_equip),
    geom = "point",
    size = 3,
    alpha = 1,
    colour = "grey30",
  ) +
  scale_shape_manual(name = NULL, values = c("mean and \n95% C.I.\nbootstrap" = 23)) + #name = NULL
  scale_fill_manual(values = c( "#e41a1c", "#377eb8"),
                    breaks = c("r5", "(land)"),
                    labels = c( "r5", "(on land)")) + 
  scale_colour_manual(values = c("#377eb8", "#e41a1c")) + #, name = NULL) +
  scale_alpha_manual(values =c(1)) +
  geom_text(data = pldatPp,
            aes(x= weave_factors(placeholder, refined, crafting_rank, Ex_rank, sailor_equip),
                y = 168.7, label = n_pp), vjust = 0.4,
            colour ="grey30",
            size = 2.5) +
#    annotation_custom(img, ymin = 166, ymax = 168, xmin = 10, xmax = 11.5) +
  guides(fill = guide_legend(order = 2, 
                             title = "Sailor Equip.", 
                             override.aes = list(size = 3, shape = 23)
                             ), 
         linetype = guide_legend(order = 3,
                                 override.aes = list(color = "#377eb8")),         
         shape = guide_legend(order = 3,
                              override.aes = list(fill = "#377eb8")),
         alpha = guide_legend(order = 5, title = NULL,
                              override.aes = list(alpha = 1, color = "grey30")),
         size = guide_legend(title = " n", order = 6),
         colour = "none") + #guide_legend(order = 5),)
  theme(plot.margin = margin(c(20, 7, 7, 7), unit = "pt"))

#interaction/summary plots
#ignorer/lump EX

#regn ut nye snit og boostrap intervals

# x = casting rank, fill/colour = sailor equip, facet refined
# y = linerange/errorbar. + lines mellom i samme farge som sailor equip

# New facet label names 
casting.labs <- c("Casting rank 3", "Casting rank 20")
names(casting.labs) <- c("r3", "r20")

ex.labs <- c("No EX", "EX rank 6")
names(ex.labs) <- c("0", "r6")

refined.labs <- c("Unrefined", "Refined")
names(refined.labs) <- c("no", "yes")

minions <- minions |>
  mutate(treatment2 = interaction(sailor_equip, refined, crafting_rank)) 

PlotMeansPp <- minions |>
  group_by(treatment2) |>
  summarise(mean_pp = mean(pp, na.rm = TRUE),
            sd_pp = sd(pp, na.rm = TRUE),
            n_pp = sum(!is.na(pp)),
            sailor_equip = unique(sailor_equip),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank))

treatments2 <- unique(PlotMeansPp$treatment2)
n_treatments2 <- length(unique(PlotMeansPp$treatment2))
### set up dataframe to collect results for plotting:
PlotdatPp <-cbind(PlotMeansPp, mean=rep(NA,n_treatments2),
                low=rep(NA, n_treatments2),
                upp=rep(NA, n_treatments2))
rownames(PlotdatPp) <-c(treatments2)


set.seed(666)
for(i in c(treatments2)){
### the boot sample and the ci's:
    boot<-boot(minions$pp[minions$treatment2==i], mean.fun, R=10000)
    ci<-boot.ci(boot,type = c("norm", "basic", "perc", 
    "bca"))
### get ci's (method: bca)
    ci[2]->PlotdatPp[i,"mean"] # mean
    ci$bca[1,4]->PlotdatPp[i,"low"]  # lower ci
    ci$bca[1,5]->PlotdatPp[i,"upp"]} # upper ci

minions |>
 ggplot(aes(x= pp, fill = sailor_equip, by = c(sailor_equip), y = after_stat(prop))) +
  geom_bar(stat = "prop", 
           position = position_dodge2(preserve = "single", padding = 0),
           linejoin = 'mitre', lineend = 'square') +
  scale_fill_manual(values = c( "#e41a1c", "#377eb8"),
                    breaks = c("r5", "(land)"),
                    labels = c( "r5", "(on land)"))  + 
  xlab("Piercing Power") +
  ylab("Frequency") +
  scale_x_continuous(breaks = c(152, 156, 160, 164, 168), minor_breaks = NULL)+
  scale_y_continuous(labels = scales::percent,
                     expand = expansion(mult = c(0, 0.05))) +
  facet_nested(vars(refined), 
               vars(crafting_rank, Ex_rank), 
               render_empty = F,
               labeller = labeller(refined = refined.labs, Ex_rank = ex.labs, crafting_rank = casting.labs)
               ) +
  guides(fill = guide_legend(title = "Sailor Equip.")) +
  theme(#strip.text.x = element_text(margin = margin(1,1,1,1, "pt")),
        #panel.spacing = unit(0.5, "lines"),
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.55, 0.8), # c(0,0) bottom left, c(1,1) top-right.
        #legend.background = element_rect(fill = "white", colour = NA)
        )
#shift_legend4(hist1)
PlotdatPp |>
  ggplot(aes(y= mean, x = crafting_rank, colour=sailor_equip)) +
    geom_line(aes(group = sailor_equip), position=position_dodge(width=.2)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=.2)) +
  facet_wrap(vars(refined), labeller = labeller(refined = refined.labs)) +
  scale_y_continuous(breaks = c(152, 156, 160, 164, 168), minor_breaks = NULL, limits = c(152, 168))+
  ylab("Piercing Power") +
  xlab("Casting rank") +
  scale_colour_manual(values = c( "#e41a1c", "#377eb8"),
                      breaks = c("r5", "(land)"),
                      labels = c( "r5", "(on land)"))  +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Sailor Equip.")) +
  theme(#strip.text.x = element_text(margin = margin(1,1,1,1, "pt")),
        #panel.spacing = unit(0.5, "lines"),
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.15, 0.8), # c(0,0) bottom left, c(1,1) top-right.
        #legend.background = element_rect(fill = "white", colour = NA)
        )

PlotdatPp |>
  ggplot(aes(y= mean, x = refined, colour=sailor_equip)) +
    geom_line(aes(group = sailor_equip), position=position_dodge(width=.2)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=.2)) +
  facet_wrap(vars(crafting_rank), labeller = labeller(crafting_rank = casting.labs)) +
  scale_y_continuous(breaks = c(152, 156, 160, 164, 168), minor_breaks = NULL, limits = c(152, 168))+
  ylab("Piercing Power") +
  xlab("Refined?") +
  scale_colour_manual(values = c( "#e41a1c", "#377eb8"),
                      breaks = c("r5", "(land)"),
                      labels = c( "r5", "(on land)")) +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Sailor Equip.")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.15, 0.8)
        )

PlotdatPp |>
  ggplot(aes(y= mean, x = sailor_equip, colour=crafting_rank)) +
    geom_line(aes(group = crafting_rank), position=position_dodge(width=.2)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=.2)) +
  facet_wrap(vars(refined), labeller = labeller(refined = refined.labs)) +
  scale_y_continuous(breaks = c(152, 156, 160, 164, 168), minor_breaks = NULL, limits = c(152, 168))+
  ylab("Piercing Power") +
  xlab("Sailor Equipment") +
  scale_colour_manual(values = c("#377eb8", "#e41a1c"), labels = c("r3", "r20")) +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Casting Rank")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.15, 0.8)
        )

PlotdatPp |>
  ggplot(aes(y= mean, x = sailor_equip, colour=refined)) +
  geom_line(aes(group = refined), position=position_dodge(width=0.2)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=0.2), 
                      alpha = 1) +
  facet_wrap(vars(crafting_rank), labeller = labeller(crafting_rank = casting.labs)) +
  scale_y_continuous(breaks = c(152, 156, 160, 164, 168), minor_breaks = NULL, limits = c(152, 168))+
  ylab("Piercing Power") +
  xlab("Sailor Equipment") +
  scale_colour_manual(values = c("#377eb8", "#e41a1c"), labels = c("unrefined", "refined")) +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Refinement")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.15, 0.8)
        )

PlotdatPp |>
  ggplot(aes(y= mean, x = crafting_rank, colour=refined)) +
  geom_line(aes(group = refined), position=position_dodge(width=0.2)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=0.2), 
                      alpha = 1) +
  facet_wrap(vars(sailor_equip)) +
  scale_y_continuous(breaks = c(152, 156, 160, 164, 168), minor_breaks = NULL, limits = c(152, 168))+
  ylab("Piercing Power") +
  xlab("Casting Rank") +
  scale_colour_manual(values = c("#377eb8", "#e41a1c"), labels = c("unrefined", "refined")) +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Refinement")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.15, 0.8), 
        )

predictionplotlmpp
anovaTable_pp

minions |>
  ggplot(aes(y = durability, x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, sailor_equip))) +
  # geom_boxplot(,
  #              outliers = F,
  #              staplewidth = 0.7,
  #              alpha = 0.5
  #              #linejoin = 'mitre', lineend = 'square'
  #              ) +
  stat_summary(
    fun = mean,
    aes(shape = "mean and \n95% C.I.\nbootstrap", fill = sailor_equip),
    geom = "point",
    size = 3,
    alpha = 1,
    color = "grey30",
  ) +
  scale_shape_manual(NULL, values = c("mean and \n95% C.I.\nbootstrap" = 23)) +
  geom_count(alpha = 0.6, shape = 16, colour = "grey70") +
  scale_size_area(breaks = c(1, 10, 25, 50, 100)) +
  ylab("Durability") +
  scale_fill_manual(values = c( "#e41a1c", "#377eb8"),
                    breaks = c("r5", "(land)"),
                    labels = c( "r5", "(on land)"))  + #na.value = "#e41a1c"
  scale_colour_manual(values = c("#377eb8", "#e41a1c"))+
  theme(axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 9, 
                                   vjust = 0.5, 
                                   margin = margin(1.5,0,1.5,0, "pt")),
        legend.key.spacing.y = unit(0.5, "pt")) + #vjust 1.5
  scale_x_discrete(guide = guide_axis_nested(
    title = NULL,
    type = "box",
    drop = FALSE
  )) +
  theme_guide(box = element_rect(colour = "grey92")) +
  annotate(
    geom = "richtext",
    size = 3.33,
    label = paste("**Factors**", "<br>refined", "<br>casting rank", "<br>EX rank", "<br>sailor equipment"),
    fill = NA,
    label.color = NA,
    x = 12.4,
    y = 105.3,
    hjust = 0,
    vjust = 1
  ) +
  coord_cartesian(expand = FALSE,
                  clip = 'off',
                  ylim = c(105, 140),
                  xlim = c(0.7, 12.3)) +
  geom_linerange(data = pldatDur, aes(
    y = mean,
    x = weave_factors(placeholder, refined, crafting_rank, Ex_rank, sailor_equip),
    ymin = low,
    ymax = upp,
    linetype = "mean and \n95% C.I.\nbootstrap", 
    colour = sailor_equip)) +
  scale_linetype_discrete(name = NULL) +
  geom_text(data = pldatDur,
            aes(x= weave_factors(placeholder, refined, crafting_rank, Ex_rank, sailor_equip),
                y = 141.5, label = n_dur), vjust = 0.4,
            colour ="grey30",
            size = 2.5) +
  guides(fill = guide_legend(order = 2, 
                             title = "Sailor Equip.", 
                             override.aes = list(size = 3, shape = 23)), 
         linetype = guide_legend(order = 3, override.aes = list(alpha = 1, color = "#377eb8")),
         shape = guide_legend(order = 3, override.aes = list(fill = "#377eb8")),
         size = guide_legend(title = " n", order = 5),
         colour = "none") +
  theme(plot.margin = margin(c(20, 7, 7, 7), unit = "pt"))
  #annotation_custom(img, ymin = 135, ymax = 140, xmin = 10, xmax = 11.5)

minions |>
 ggplot(aes(x = durability, fill = sailor_equip, by = c(sailor_equip))) + #y = after_stat(prop)
  geom_histogram(aes(x= durability, after_stat(density*width)), binwidth = 8, #stat = "prop", 
           position = position_dodge2(preserve = "single", padding = 0),
           linejoin = 'mitre', lineend = 'square') +
  scale_fill_manual(values = c( "#e41a1c", "#377eb8"),
                    breaks = c("r5", "(land)"),
                    labels = c( "r5", "(on land)"))  + 
  xlab("Durability") +
  ylab("Frequency") +
  #scale_x_continuous(breaks = c(152, 156, 160, 164, 168))+
  scale_y_continuous(labels = scales::percent,
                         expand = expansion(mult = c(0, 0.05))) +
  facet_nested(vars(refined), 
               vars(crafting_rank, Ex_rank), 
               render_empty = F,
               labeller = labeller(refined = refined.labs, Ex_rank = ex.labs, crafting_rank = casting.labs)
               ) +
  guides(fill = guide_legend(title = "Sailor Equip.")) +
  theme(#strip.text.x = element_text(margin = margin(1,1,1,1, "pt")),
        #panel.spacing = unit(0.5, "lines"),
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.55, 0.8) # c(0,0) bottom left, c(1,1) top-right.
        #legend.background = element_rect(fill = "white", colour = NA)
        ) 

set.seed(666)
PlotMeansDur <- minions |>
  group_by(treatment2) |>
  summarise(mean_dur = mean(durability, na.rm = TRUE),
            sd_dur = sd(durability, na.rm = TRUE),
            n_dur =  sum(!is.na(durability)),
            sailor_equip = unique(sailor_equip),
            refined = unique(refined),
            crafting_rank = unique(crafting_rank))

PlotdatDur <-cbind(PlotMeansDur, mean=rep(NA,n_treatments2),
                 low=rep(NA, n_treatments2),
                 upp=rep(NA, n_treatments2))

rownames(PlotdatDur) <-c(treatments2)


for(i in c(treatments2)){
    boot<-boot(minions$durability[minions$treatment2 == i], mean.fun, R = 10000)
    ci<-boot.ci(boot,type = c("norm", "basic", "perc", 
    "bca"))
    ci[2] -> PlotdatDur[i, "mean"] # mean
    ci$bca[1, 4]->PlotdatDur[i, "low"]  # lower ci
    ci$bca[1, 5]->PlotdatDur[i, "upp"]
    } # upper ci


PlotdatDur |>
  ggplot(aes(y= mean, x = crafting_rank, colour=sailor_equip)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=0)) +
  geom_line(aes(group = sailor_equip), position=position_dodge(width=0)) +
  facet_wrap(vars(refined), labeller = labeller(refined = refined.labs)) +
  scale_y_continuous(minor_breaks = NULL, limits = c(105, 140))+
  ylab("Durability") +
  xlab("Casting rank") +
  scale_colour_manual(values = c( "#e41a1c", "#377eb8"),
                      breaks = c("r5", "(land)"),
                      labels = c( "r5", "(on land)"))  +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Sailor Equip.")) +
  theme(#strip.text.x = element_text(margin = margin(1,1,1,1, "pt")),
        #panel.spacing = unit(0.5, "lines"),
        #panel.grid.major.x = element_blank(),
        panel.border = element_rect(color = "white", fill = NA), 
        strip.clip = "off",
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.65, 0.3), # c(0,0) bottom left, c(1,1) top-right.
        #legend.background = element_rect(fill = "white", colour = NA)
        ) #+
 # panel_border(colour = "white")

PlotdatDur |>
  ggplot(aes(y= mean, x = refined, colour=sailor_equip)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=0)) +
  geom_line(aes(group = sailor_equip), position=position_dodge(width=0)) +
  facet_wrap(vars(crafting_rank), labeller = labeller(crafting_rank = casting.labs)) +
  scale_y_continuous(minor_breaks = NULL, limits = c(105, 140))+
  ylab("Durability") +
  xlab("Refined?") +
  scale_colour_manual(values = c( "#e41a1c", "#377eb8"),
                      breaks = c("r5", "(land)"),
                      labels = c( "r5", "(on land)"))  +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Sailor Equip.")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.65, 0.3)
        )

PlotdatDur |>
  ggplot(aes(y= mean, x = sailor_equip, colour=crafting_rank)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=0)) +
  geom_line(aes(group = crafting_rank), position=position_dodge(width=0)) +
  facet_wrap(vars(refined), labeller = labeller(refined = refined.labs)) +
  scale_y_continuous(minor_breaks = NULL, limits = c(105, 140))+
  ylab("Durability") +
  xlab("Sailor Equipment") +
  scale_colour_manual(values = c("#377eb8", "#e41a1c"), labels = c("r3", "r20")) +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Casting rank")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.65, 0.3)
        )

PlotdatDur |>
  ggplot(aes(y= mean, x = sailor_equip, colour=refined)) +
  geom_line(aes(group = refined), position=position_dodge(width=0.2)) +
  geom_pointrange(aes(ymin= low, ymax= upp,
                      linetype = "mean and\n95% C.I.\nbootstrap"), 
                      #linewidth = 1, 
                      position=position_dodge(width=0.2), 
                      alpha = 1) +
  facet_wrap(vars(crafting_rank), labeller = labeller(crafting_rank = casting.labs)) +
  scale_y_continuous(minor_breaks = NULL, limits = c(105, 140))+
  ylab("Durability") +
  xlab("Sailor Equipment") +
  scale_colour_manual(values = c("#377eb8", "#e41a1c"), labels = c("unrefined", "refined")) +
  scale_linetype_discrete(name = NULL) +
  guides(colour = guide_legend(title = "Refinement")) +
  theme(
        panel.border = element_rect(color = "white", fill = NA), 
        strip.background = element_rect(color = "white"),
        legend.position = "inside",
        legend.position.inside = c(0.65, 0.3)
        )

#lag et plot med EX også for å vise at det er ingen forskjell (bruk colour)
#må lage nye bootstrap CI da

anovaTable_dura
cor.test(minions$pp, minions$durability, method = "spearman")

minions |>
  ggplot(aes(x = pp, y = durability)) +
  geom_count() +
  scale_size_area() 

minions |>
  split(minions$treatment) |>
  map(function(df) cor.test(df$pp, df$durability, method = "spearman")) |>
  map_dbl("estimate")

minions |> 
  filter(treatment == "J") |>
  ggplot(aes(x = pp, y = durability)) +
  geom_count() +
  scale_size_area() 

predictionplotlm
modelSummaryGLM

AICglm
predictionplotGlm

modelSummaryPp

AICpp
predictionplotlmpp
modelSummaryDura

AICdura
predictionplotlm

#PERMUTATION TESTS OF BATCH EFFECTS
# sjekk om variasjonen observert i batches er mindre enn tilfeldige batches av 20 trukket fra data (innen treatment) og variasjonen mellom batcher større.

#må holde orden på sample size i hver batch, og trekke tilsvarende antall tilfedlig

#variance per batch
obsVariances<- minions |>
  group_by(treatment, batch_id) |>
  dplyr::summarise(variancePp = var(pp, na.rm = T),
                   sdPp = sd(pp, na.rm = T),
             varianceDur = var(durability, na.rm = T),
             sdDur = sd(durability, na.rm = T),
             valid_n = sum(!is.na(pp)))
# mean(obsVariances$varianceDur)
# hist(obsVariances$varianceDur)
# mean(obsVariances$variancePp)
# hist(obsVariances$variancePp)
# mean(obsVariances$sdDur)
# hist(obsVariances$sdDur)
# mean(obsVariances$sdPp)
# hist(obsVariances$sdPp)
# 
# mean(minions$sd_batchDura, na.rm =T)
# mean(minions$sd_batchPp, na.rm =T)

#group by treatment, slice sample pp but keep the other variables
#size needs to be a named vector
simVariances <- replicate(10000, simplify = F,
                          Rsampling::within_columns(reducedMinions, 
                                          cols = c("durability", "pp"), stratum = reducedMinions$treatment) |>
  group_by(treatment, batch_id) |>
  dplyr::summarise(variancePp = var(pp, na.rm = T),
             varianceDur = var(durability, na.rm = T),
             sdPp = sd(pp, na.rm = T),
             sdDur = sd(durability, na.rm = T),
             valid_n = sum(!is.na(pp))) 
  |> ungroup() |>
    dplyr::summarise(meanVarPp = mean(variancePp),
              meanVarDur = mean(varianceDur),
              meanSdPp = mean(sdPp),
              meanSdDur = mean(sdDur))
  )

simVariances <- as.data.frame(do.call(rbind, simVariances ))

save(simVariances, file = "simVariances.RData")

load(file = "simVariances.RData")

p_valuePpSd <- 
  round(length(which(simVariances$meanSdPp >= mean(obsVariances$sdPp)))/length(simVariances[,1])*2, 2)

simVariances |>
  ggplot(aes(x=meanSdPp)) +
  geom_histogram(aes(fill = "simulated"), bins = 30) +
  geom_vline(aes(linetype = "observed", xintercept =  mean(obsVariances$sdPp)), colour = "#377eb8", linewidth = 1.5) +
  guides(linetype = guide_legend(title = "SD", order = 1),
         fill = guide_legend(title = NULL, order = 2)) +
    scale_fill_manual(values = c("#A6CEE3")) +
  xlab("Mean Standard Deviation (SD) of Piercing Power \namong batches (permuted within treatments)") +
  ylab("Frequency") +
    annotate(geom = "text",
           label = paste0("n = ", length(reducedMinions$pp) , " minion(4) cannons",
                          "\nin ",length(obsVariances$sdPp), " batches",
                          "\npermutations within treatments",
                          "\n10,000 simulations",
                          "\np = ", p_valuePpSd),
           x = 5.65, 
           y = 850,
           hjust = 0,
           size = 3
    ) +
    scale_y_continuous(
    expand = expansion(mult = c(0, 0.05))) +
  theme(
        legend.position= "inside",
        legend.position.inside = c(1/7, 1/2),
        legend.margin = NULL)+
  coord_cartesian(clip = 'off')
p_valueDuraSd <- 
  round(length(which(simVariances$meanSdDur >= mean(obsVariances$sdDur)))/length(simVariances[,1])*2, 2)

simVariances |>
  ggplot(aes(x=meanSdDur)) +
  geom_histogram(aes(fill = "simulated"), bins = 30) +
  geom_vline(aes(linetype = "observed", xintercept =  mean(obsVariances$sdDur)), colour = "#377eb8", linewidth = 1.5) +
  guides(linetype = guide_legend(title = "SD", order = 1),
         fill = guide_legend(title = NULL, order = 2)) +
    scale_fill_manual(values = c("#A6CEE3")) +
  xlab("Mean Standard Deviation (SD) of Durability \namong batches (permuted within treatments)") +
  ylab("Frequency") +
    annotate(geom = "text",
           label = paste0("n = ", length(reducedMinions$pp) , " minion(4) cannons",
                          "\nin ",length(obsVariances$sdPp), " batches",
                          "\npermutations within treatments",
                          "\n10,000 simulations",
                          "\np = ", p_valueDuraSd),
           x = 7.25, 
           y = 850,
           hjust = 0,
           size = 3
    ) +
  scale_y_continuous(
    expand = expansion(mult = c(0, 0.05))
  ) +
  theme(
        legend.position= "inside",
        legend.position.inside = c(3/4, 0.7),
        legend.margin = NULL)+
  coord_cartesian(clip = 'off')
  
pkg_sesh <- sessioninfo::session_info()
quarto_version <- system("quarto --version", intern = TRUE)
pkg_sesh$platform$quarto <- paste(
  system("quarto --version", intern = TRUE)
  )
pkg_sesh
library(downloadthis)

minions |>
  download_this(
    output_name = "minionData",
    output_extension = ".csv",
    button_label = "Download Minion data as csv",
    button_type = "primary",
    has_icon = TRUE,
    icon = "fa fa-save",
    class = "hvr-sweep-to-left"
  )
additional |>
  download_this(
    output_name = "additionalMinionData",
    output_extension = ".csv",
    button_label = "Download additional Minion data as csv",
    button_type = "primary",
    has_icon = TRUE,
    icon = "fa fa-save",
    class = "hvr-sweep-to-left"
  )

Software information

Session Info

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31 ucrt)
 os       Windows 10 x64 (build 19045)
 system   x86_64, mingw32
 ui       RTerm
 language (EN)
 collate  nb_NO.utf8
 ctype    nb_NO.utf8
 tz       Europe/Oslo
 date     2025-02-25
 pandoc   3.2 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
 quarto   1.6.40

─ Packages ───────────────────────────────────────────────────────────────────
 package           * version    date (UTC) lib source
 abind               1.4-8      2024-09-12 [1] CRAN (R 4.4.1)
 afex              * 1.4-1      2024-09-01 [1] CRAN (R 4.4.1)
 AICcmodavg          2.3-3      2023-11-16 [1] CRAN (R 4.3.1)
 arm                 1.14-4     2024-04-01 [1] CRAN (R 4.3.3)
 backports           1.5.0      2024-05-23 [1] CRAN (R 4.3.3)
 base64enc           0.1-3      2015-07-28 [1] CRAN (R 4.3.0)
 bayestestR          0.15.2     2025-02-07 [1] CRAN (R 4.4.2)
 binom               1.1-1.1    2022-05-02 [1] CRAN (R 4.3.0)
 bitops              1.0-9      2024-10-03 [1] CRAN (R 4.4.1)
 boot              * 1.3-31     2024-08-28 [1] CRAN (R 4.4.2)
 boot.pval           0.6        2025-01-14 [1] CRAN (R 4.4.2)
 broom               1.0.7      2024-09-26 [1] CRAN (R 4.4.1)
 broom.helpers       1.19.0     2025-01-29 [1] CRAN (R 4.4.2)
 broom.mixed         0.2.9.6    2024-10-15 [1] CRAN (R 4.4.1)
 car                 3.1-3      2024-09-27 [1] CRAN (R 4.4.1)
 carData             3.0-5      2022-01-06 [1] CRAN (R 4.3.0)
 cards               0.5.0      2025-02-17 [1] CRAN (R 4.4.2)
 cardx               0.2.3      2025-02-18 [1] CRAN (R 4.4.2)
 caTools             1.18.3     2024-09-04 [1] CRAN (R 4.4.1)
 cellranger          1.1.0      2016-07-27 [1] CRAN (R 4.3.0)
 checkmate           2.3.2      2024-07-29 [1] CRAN (R 4.4.1)
 class               7.3-23     2025-01-01 [1] CRAN (R 4.4.2)
 cli                 3.6.4      2025-02-13 [1] CRAN (R 4.4.2)
 coda                0.19-4.1   2024-01-31 [1] CRAN (R 4.3.2)
 codetools           0.2-20     2024-03-31 [2] CRAN (R 4.4.2)
 colorspace          2.1-1      2024-07-26 [1] CRAN (R 4.4.1)
 commonmark          1.9.2      2024-10-04 [1] CRAN (R 4.4.1)
 cowplot           * 1.1.3      2024-01-22 [1] CRAN (R 4.3.2)
 data.table          1.16.4     2024-12-06 [1] CRAN (R 4.4.2)
 datawizard          1.0.0      2025-01-10 [1] CRAN (R 4.4.2)
 DEoptimR            1.1-3-1    2024-11-23 [1] CRAN (R 4.4.1)
 digest              0.6.37     2024-08-19 [1] CRAN (R 4.4.1)
 distributional      0.5.0      2024-09-17 [1] CRAN (R 4.4.1)
 doParallel          1.0.17     2022-02-07 [1] CRAN (R 4.3.0)
 dplyr             * 1.1.4      2023-11-17 [1] CRAN (R 4.3.2)
 e1071               1.7-16     2024-09-16 [1] CRAN (R 4.4.1)
 effectsize          1.0.0      2024-12-10 [1] CRAN (R 4.4.1)
 emmeans             1.10.7     2025-01-31 [1] CRAN (R 4.4.2)
 estimability        1.5.1      2024-05-12 [1] CRAN (R 4.3.3)
 evaluate            1.0.3      2025-01-10 [1] CRAN (R 4.4.2)
 farver              2.1.2      2024-05-13 [1] CRAN (R 4.3.3)
 fastmap             1.2.0      2024-05-15 [1] CRAN (R 4.3.3)
 forcats             1.0.0      2023-01-29 [1] CRAN (R 4.3.0)
 foreach             1.5.2      2022-02-02 [1] CRAN (R 4.3.0)
 Formula             1.2-5      2023-02-24 [1] CRAN (R 4.3.0)
 furrr               0.3.1      2022-08-15 [1] CRAN (R 4.3.0)
 future              1.34.0     2024-07-29 [1] CRAN (R 4.4.1)
 future.apply        1.11.3     2024-10-27 [1] CRAN (R 4.4.1)
 generics            0.1.3      2022-07-05 [1] CRAN (R 4.3.0)
 ggdist            * 3.3.2      2024-03-05 [1] CRAN (R 4.3.3)
 ggh4x             * 0.3.0      2024-12-15 [1] CRAN (R 4.4.1)
 ggplot2           * 3.5.1      2024-04-23 [1] CRAN (R 4.3.3)
 ggResidpanel        0.3.0      2019-05-31 [1] CRAN (R 4.4.2)
 ggstats           * 0.8.0      2025-01-07 [1] CRAN (R 4.4.1)
 ggtext            * 0.1.2      2022-09-16 [1] CRAN (R 4.3.0)
 globals             0.16.3     2024-03-08 [1] CRAN (R 4.3.3)
 glue                1.8.0      2024-09-30 [1] CRAN (R 4.4.1)
 gridtext            0.1.5      2022-09-16 [1] CRAN (R 4.3.0)
 gt                * 0.11.1     2024-10-04 [1] CRAN (R 4.4.1)
 gtable              0.3.6      2024-10-25 [1] CRAN (R 4.4.1)
 gtsummary         * 2.1.0      2025-02-19 [1] CRAN (R 4.4.2)
 haven               2.5.4      2023-11-30 [1] CRAN (R 4.3.2)
 hms                 1.1.3      2023-03-21 [1] CRAN (R 4.3.0)
 htmltools           0.5.8.1    2024-04-04 [1] CRAN (R 4.3.3)
 htmlwidgets         1.6.4      2023-12-06 [1] CRAN (R 4.3.2)
 httr                1.4.7      2023-08-15 [1] CRAN (R 4.3.1)
 insight             1.0.2      2025-02-06 [1] CRAN (R 4.4.2)
 iterators           1.0.14     2022-02-05 [1] CRAN (R 4.3.0)
 jsonlite            1.9.0      2025-02-19 [1] CRAN (R 4.4.2)
 jtools              2.3.0      2024-08-25 [1] CRAN (R 4.4.2)
 knitr               1.49       2024-11-08 [1] CRAN (R 4.4.1)
 labeling            0.4.3      2023-08-29 [1] CRAN (R 4.3.1)
 labelled            2.14.0     2025-01-08 [1] CRAN (R 4.4.1)
 lattice             0.22-6     2024-03-20 [1] CRAN (R 4.3.3)
 lazyeval            0.2.2      2019-03-15 [1] CRAN (R 4.3.0)
 legendry          * 0.2.0.9000 2025-02-21 [1] Github (teunbrand/legendry@83ee986)
 lifecycle           1.0.4      2023-11-07 [1] CRAN (R 4.3.1)
 listenv             0.9.1      2024-01-29 [1] CRAN (R 4.3.2)
 lme4              * 1.1-36     2025-01-11 [1] CRAN (R 4.4.2)
 lmerTest            3.1-3      2020-10-23 [1] CRAN (R 4.3.0)
 lmtest              0.9-40     2022-03-21 [1] CRAN (R 4.3.0)
 magick              2.8.5      2024-09-20 [1] CRAN (R 4.4.1)
 magrittr            2.0.3      2022-03-30 [1] CRAN (R 4.3.0)
 marginaleffects     0.25.0     2025-02-01 [1] CRAN (R 4.4.2)
 markdown            1.13       2024-06-04 [1] CRAN (R 4.3.3)
 MASS                7.3-64     2025-01-04 [1] CRAN (R 4.4.2)
 Matrix            * 1.7-2      2025-01-23 [1] CRAN (R 4.4.2)
 minqa               1.2.8      2024-08-17 [1] CRAN (R 4.4.1)
 minty               0.0.5      2025-01-07 [1] CRAN (R 4.4.2)
 modelsummary      * 2.3.0      2025-02-02 [1] CRAN (R 4.4.2)
 multcomp            1.4-28     2025-01-29 [1] CRAN (R 4.4.2)
 munsell             0.5.1      2024-04-01 [1] CRAN (R 4.3.3)
 mvtnorm             1.3-3      2025-01-10 [1] CRAN (R 4.4.2)
 nlme                3.1-167    2025-01-27 [1] CRAN (R 4.4.2)
 nloptr              2.1.1      2024-06-25 [1] CRAN (R 4.3.3)
 numDeriv            2016.8-1.1 2019-06-06 [1] CRAN (R 4.3.0)
 opdisDownsampling   1.0.1      2024-04-15 [1] CRAN (R 4.4.2)
 pander              0.6.5      2022-03-18 [1] CRAN (R 4.3.0)
 parallelly          1.42.0     2025-01-30 [1] CRAN (R 4.4.2)
 parameters          0.24.1     2025-01-14 [1] CRAN (R 4.4.2)
 PASWR2              1.0.5      2021-09-04 [1] CRAN (R 4.4.2)
 patchwork         * 1.3.0      2024-09-16 [1] CRAN (R 4.4.1)
 pbmcapply           1.5.1      2022-04-28 [1] CRAN (R 4.4.0)
 performance         0.13.0     2025-01-15 [1] CRAN (R 4.4.2)
 pillar              1.10.1     2025-01-07 [1] CRAN (R 4.4.1)
 pkgconfig           2.0.3      2019-09-22 [1] CRAN (R 4.3.0)
 plotly              4.10.4     2024-01-13 [1] CRAN (R 4.3.2)
 plyr                1.8.9      2023-10-02 [1] CRAN (R 4.3.1)
 pracma              2.4.4      2023-11-10 [1] CRAN (R 4.3.2)
 proxy               0.4-27     2022-06-09 [1] CRAN (R 4.3.0)
 purrr             * 1.0.4      2025-02-05 [1] CRAN (R 4.4.2)
 qqconf              1.3.2      2023-04-14 [1] CRAN (R 4.3.0)
 qqplotr             0.0.6      2023-01-25 [1] CRAN (R 4.4.2)
 quadprog            1.5-8      2019-11-20 [1] CRAN (R 4.3.0)
 R6                  2.6.1      2025-02-15 [1] CRAN (R 4.4.2)
 ragg                1.3.3      2024-09-11 [1] CRAN (R 4.4.1)
 rbibutils           2.3        2024-10-04 [1] CRAN (R 4.4.1)
 RColorBrewer        1.1-3      2022-04-03 [1] CRAN (R 4.3.0)
 Rcpp                1.0.14     2025-01-12 [1] CRAN (R 4.4.2)
 Rdpack              2.6.2      2024-11-15 [1] CRAN (R 4.4.2)
 readODS             2.3.2      2025-01-13 [1] CRAN (R 4.4.2)
 reformulas          0.4.0      2024-11-03 [1] CRAN (R 4.4.1)
 reshape2            1.4.4      2020-04-09 [1] CRAN (R 4.3.0)
 rlang               1.1.5      2025-01-17 [1] CRAN (R 4.4.2)
 rmarkdown           2.29       2024-11-04 [1] CRAN (R 4.4.1)
 robustbase          0.99-4-1   2024-09-27 [1] CRAN (R 4.4.1)
 rstudioapi          0.17.1     2024-10-22 [1] CRAN (R 4.4.1)
 sandwich          * 3.1-1      2024-09-15 [1] CRAN (R 4.4.1)
 sass                0.4.9      2024-03-15 [1] CRAN (R 4.3.3)
 scales              1.3.0      2023-11-28 [1] CRAN (R 4.3.2)
 sessioninfo         1.2.3      2025-02-05 [1] CRAN (R 4.4.2)
 stringi             1.8.4      2024-05-06 [1] CRAN (R 4.3.3)
 stringr             1.5.1      2023-11-14 [1] CRAN (R 4.3.1)
 survival            3.8-3      2024-12-17 [1] CRAN (R 4.4.2)
 systemfonts         1.2.1      2025-01-20 [1] CRAN (R 4.4.2)
 tables              0.9.31     2024-08-29 [1] CRAN (R 4.4.1)
 textshaping         1.0.0      2025-01-20 [1] CRAN (R 4.4.2)
 TH.data             1.1-3      2025-01-17 [1] CRAN (R 4.4.2)
 tibble              3.2.1      2023-03-20 [1] CRAN (R 4.3.0)
 tidyr             * 1.3.1      2024-01-24 [1] CRAN (R 4.3.2)
 tidyselect          1.2.1      2024-03-11 [1] CRAN (R 4.3.3)
 twosamples          2.0.1      2023-06-23 [1] CRAN (R 4.4.2)
 tzdb                0.4.0      2023-05-12 [1] CRAN (R 4.3.0)
 unmarked            1.5.0      2025-02-10 [1] CRAN (R 4.4.2)
 vctrs               0.6.5      2023-12-01 [1] CRAN (R 4.3.2)
 VGAM                1.1-13     2025-02-12 [1] CRAN (R 4.4.2)
 viridisLite         0.4.2      2023-05-02 [1] CRAN (R 4.3.0)
 withr               3.0.2      2024-10-28 [1] CRAN (R 4.4.1)
 xfun                0.51       2025-02-19 [1] CRAN (R 4.4.2)
 xml2                1.3.6      2023-12-04 [1] CRAN (R 4.3.2)
 xtable              1.8-4      2019-04-21 [1] CRAN (R 4.3.0)
 yaml                2.3.10     2024-07-26 [1] CRAN (R 4.4.1)
 zip                 2.3.2      2025-02-01 [1] CRAN (R 4.4.2)
 zoo                 1.8-12     2023-04-13 [1] CRAN (R 4.4.2)

 [1] C:/R/RLibs
 [2] C:/Program Files/R/R-4.4.2/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────

Raw data