set.seed(123456789)
library(MASS)
library(readxl)
library(EnvStats)
##
## Attaching package: 'EnvStats'
## The following object is masked from 'package:MASS':
##
## boxcox
## The following objects are masked from 'package:stats':
##
## predict, predict.lm
## The following object is masked from 'package:base':
##
## print.default
presences.num <- (read.csv("presences.num.csv"))[ , -1] #Remove column 1, row numbers
presences_lessLL <- presences.num[ , -23] # Remove Lucy, now column 23, named "L3"
presences_lessLL_T <- t(presences_lessLL)
dim(presences_lessLL)
## [1] 25 37
colnames(presences_lessLL_T) <- c(" 1", " 2", " 3", " 4", " 5", " 6", " 7", " 8", " 9", 10:25)
presences_lessLL_T
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## A1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## A2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
## A3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0
## A4 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 1 1 0 0 0 0 0 0
## B1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0
## B2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## C1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## C2 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0
## C3 0 0 0 1 1 0 0 1 0 0 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0
## C4 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## C5 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 1 0 1 1 1 0
## E1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## E2 1 0 1 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## E3 0 0 1 0 0 0 0 1 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 0 1
## J1 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 1 0
## J2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0
## J3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## J4 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0
## K1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## K2 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## L1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0
## L2 1 1 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## M1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 1 1
## M2 1 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 1 0 0
## M3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
## N1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1
## P1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0 1 0 0 1 1 0 0
## R1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0
## S1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0
## S2 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
## S3 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0
## S4 1 1 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## V1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0
## V2 0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0
## V3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0
## Y1 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 0 0 0 1 1 1 0
## Y2 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0
sum(presences_lessLL) / (25 * 37)
## [1] 0.1556757
nurse_totals <- colSums(presences_lessLL)
nurse_totals[1 : 22]
## A1 A2 A3 A4 B1 B2 C1 C2 C3 C4 C5 E1 E2 E3 J1 J2 J3 J4 K1 K2 L1 L2
## 1 1 2 5 5 2 2 3 7 2 7 3 5 7 6 2 1 3 3 2 4 6
nurse_totals[23 : 37]
## M1 M2 M3 N1 P1 R1 S1 S2 S3 S4 V1 V2 V3 Y1 Y2
## 7 7 2 3 6 3 2 3 4 4 5 6 2 6 5
ebeta.fit <- ebeta(nurse_totals/25, method = "mle")
ebeta.fit
## $distribution
## [1] "Beta"
##
## $sample.size
## [1] 37
##
## $parameters
## shape1 shape2
## 3.084949 16.773344
##
## $n.param.est
## [1] 2
##
## $method
## [1] "mle"
##
## $data.name
## [1] "nurse_totals/25"
##
## $bad.obs
## [1] 0
##
## attr(,"class")
## [1] "estimate"
alpha <- ebeta.fit$parameters["shape1"]
beta <- ebeta.fit$parameters["shape2"]
rbetas <- rbeta(25, alpha, beta)
truehist(rbetas, h = 0.05, xlim = c(0, 1), ylim = c(0, 7), prob = TRUE,
main = "Histogram of a sample of size 25 \n from another fitted beta distribution",
sub = "alpha = 3.08, beta = 16.77",
xlab = "Nurse specific probability of shift rostering", ylab = "Probability density")
p <- (0:1000)/1000
lines(p, dbeta(p, alpha, beta), col = "red", lwd = 2)
lines(x = c(0, 1), y = c(0, 0)); lines(x = c(0, 0), y = c(0, 7))
