### R code from vignette source 'ReproducingN2025.Rnw'

###################################################
### code chunk number 1: ReproducingN2025.Rnw:76-82
###################################################
# Call library and data. Give short names to data objects.
library(apc)
data <- data.Swiss.suicides()
data.list.f<-data$data.list.f
data.list.m<-data$data.list.m
F_index<-data$v.F_index


###################################################
### code chunk number 2: ReproducingN2025.Rnw:87-90
###################################################
# Show some of the data 
objects(data)
(data.list.f)$response[1:5,1:5]


###################################################
### code chunk number 3: ReproducingN2025.Rnw:97-99 (eval = FALSE)
###################################################
## # Various plots of data sums
## apc.plot.data.sums(data.list.f,"c",scale.rate=100000)


###################################################
### code chunk number 4: ReproducingN2025.Rnw:103-117
###################################################
# Function for plotting two variables in one plot
plot.both<- function(y.f,y.m,main,position){
x<-1:length(y.f)
labels<-names(y.f)
xlim<-c(1,length(x))
ylim<-c(0,max(y.f,y.m))
plot(NULL,xaxt="n",xlab="",ylab="",xlim=xlim,ylim=ylim,main=main)
axis(side=1,at=x,labels)
legend(position,col=c("black","blue"),lty=c(1,1),pch=c(19,1),legend=c("women","men"))
 lines(x,y.f,col="black")
points(x,y.f,col="black",pch=19) 
 lines(x,y.m,col="blue" )         
points(x,y.m,col="blue" ,pch=1 )
}


###################################################
### code chunk number 5: ReproducingN2025.Rnw:120-125
###################################################
# Compute crude rates
crude.age.f<-100000*rowSums(data.list.f$response)/rowSums(data.list.f$dose)
crude.age.m<-100000*rowSums(data.list.m$response)/rowSums(data.list.m$dose)
crude.per.f<-100000*colSums(data.list.f$response)/colSums(data.list.f$dose)
crude.per.m<-100000*colSums(data.list.m$response)/colSums(data.list.m$dose)


###################################################
### code chunk number 6: ReproducingN2025.Rnw:128-130
###################################################
# Plot of data sums by age
plot.both(crude.age.f,crude.age.m,"Crude rate by age","topleft")


###################################################
### code chunk number 7: ReproducingN2025.Rnw:132-134
###################################################
# Plot of data sums by period
plot.both(crude.per.f,crude.per.m,"Crude rate by period","topright")


###################################################
### code chunk number 8: ReproducingN2025.Rnw:141-155
###################################################
# Analysis of variance, women
apc.table.f<-apc.fit.table(data$data.list.f,"log.normal.rates")
Table2<-apc.table.f[1:4,c(1,2,6,4,7)]
# Add scale estimates
sigma.f<-sqrt(c(apc.fit.model(data$data.list.f,"log.normal.rates","APC")$s2,
apc.fit.model(data$data.list.f,"log.normal.rates","AP" )$s2,
apc.fit.model(data$data.list.f,"log.normal.rates","AC" )$s2,
apc.fit.model(data$data.list.f,"log.normal.rates","PC" )$s2))
sigma<-sigma.f
Table2<-cbind(Table2,sigma)
# Increase precision on p.values
Table2[,5]<-pf(as.numeric(Table2[,3]),as.numeric(Table2[,4]),as.numeric(Table2[,2]),lower.tail=FALSE)
print("Table 2: Analysis of variance, women")
Table2


###################################################
### code chunk number 9: ReproducingN2025.Rnw:158-162
###################################################
# Normality test, women			   
fit.apc.f	<- apc.fit.model(data$data.list.f,"log.normal.rates","APC")
print("Normality test, women")
apc.test.normal.residuals(fit.apc.f,remove.zeros=TRUE)[c(1,6),]


###################################################
### code chunk number 10: ReproducingN2025.Rnw:169-183
###################################################
# Analysis of variance, men
apc.table.m<-apc.fit.table(data$data.list.m,"log.normal.rates")
Table2<-apc.table.m[1:4,c(1,2,6,4,7)]
# Add scale estimates
sigma.m<-sqrt(c(apc.fit.model(data$data.list.m,"log.normal.rates","APC")$s2,
apc.fit.model(data$data.list.m,"log.normal.rates","AP" )$s2,
apc.fit.model(data$data.list.m,"log.normal.rates","AC" )$s2,
apc.fit.model(data$data.list.m,"log.normal.rates","PC" )$s2))
sigma<-sigma.m
Table2<-cbind(Table2,sigma)
# Increase precision on p.values
Table2[,5]<-pf(as.numeric(Table2[,3]),as.numeric(Table2[,4]),as.numeric(Table2[,2]),lower.tail=FALSE)
print("Table 2: Analysis of variance, men")
Table2


###################################################
### code chunk number 11: ReproducingN2025.Rnw:186-190
###################################################
# Normality test, men
fit.apc.m	<- apc.fit.model(data$data.list.m,"log.normal.rates","APC")
print("Normality test, men")
apc.test.normal.residuals(fit.apc.m,remove.zeros=TRUE)[c(1,6),]


###################################################
### code chunk number 12: ReproducingN2025.Rnw:198-206 (eval = FALSE)
###################################################
## # Analysis of variance, both
## table.ols<-apc.fit.table.2s(data.list.f,data.list.m,"log.normal.rates","difference","APC","APC")
## s.f	<- sqrt(apc.fit.model(data.list.f,"log.normal.rates","APC")$s2)
## s.m	<- sqrt(apc.fit.model(data.list.m,"log.normal.rates","APC")$s2)
## table.gls	<- apc.fit.table.2s(data.list.f,data.list.m,"gls.log.normal.rates","difference","APC","APC",c(s.f/s.m,1))
## Table4 <- cbind(table.ols[1:4,c(1,2,4,6,7,10)],table.gls[1:4,10])
## colnames(Table4)[c(1,6:7)] <- c("-2logL_OLS","sigma_OLS","sigma_GLS")
## Table4


###################################################
### code chunk number 13: ReproducingN2025.Rnw:210-217 (eval = FALSE)
###################################################
## # Normality tests, both
## fit.apc.OLS	<- apc.fit.model.2s(data.list.f,data.list.m,	"log.normal.rates","APC","APC")
## fit.apc.GLS	<- apc.fit.model.2s(data.list.f,data.list.m,"gls.log.normal.rates","APC","APC",c(s.f/s.m,1))
## print("Normality test, two-sample, OLS")
## apc.test.normal.residuals(fit.apc.OLS,remove.zeros=TRUE)[c(1,6),]
## print("Normality test, two-sample, GLS")
## apc.test.normal.residuals(fit.apc.GLS,remove.zeros=TRUE)[c(1,6),]


###################################################
### code chunk number 14: ReproducingN2025.Rnw:260-268
###################################################
# Plot 1-sample fits jointly. Macro effects.
par(mfrow=c(3,2),oma=c(0,0,2,0),mar=c(4,2,2,2)+0.1)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="macro",which.plot=1)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="macro",which.plot=4)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="macro",which.plot=2)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="macro",which.plot=5)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="macro",which.plot=3)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="macro",which.plot=6)


###################################################
### code chunk number 15: ReproducingN2025.Rnw:279-285
###################################################
# Plot 1-sample fits jointly. Micro effects.
par(mfrow=c(2,2),oma=c(0,0,2,0),mar=c(4,2,2,2)+0.1)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="micro",which.plot=2)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="micro",which.plot=5)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="micro",which.plot=3)
apc.plot.fit.2s(fit.apc.f,fit.apc.m,type="micro",which.plot=6)


###################################################
### code chunk number 16: ReproducingN2025.Rnw:294-297
###################################################
# Detrending F_index
x.F	<- 1950:2007
F_index.detrend	<- F_index-F_index[1] - (0:57)/57*(F_index[58]-F_index[1])


###################################################
### code chunk number 17: ReproducingN2025.Rnw:301-306
###################################################
# Plot F_index
plot(x.F,F_index,type="l",ylim=c(-0.2,0.8),xlab="",ylab="")
lines(x.F,-F_index.detrend,lty=2)
title(main="F-index")
legend("topright",legend=c("Raw","Detrended, *(-1)"),lty=c(1,2))	


###################################################
### code chunk number 18: ReproducingN2025.Rnw:314-321 (eval = FALSE)
###################################################
## # Analysis of variance, period replaced by time series
## s.f	<- sqrt(apc.fit.model(data.list.f,"log.normal.rates","APC")$s2)
## s.m	<- sqrt(apc.fit.model(data.list.m,"log.normal.rates","APC")$s2)
## table.gls.ts.apc<-apc.fit.table.2s(data.list.f,data.list.m,"gls.log.normal.rates","difference","APC","APC",c(s.f/s.m,1),time.series=F_index)
## table.ols.ts.apc<-apc.fit.table.2s(data.list.f,data.list.m,"log.normal.rates","difference","APC","APC",time.series=F_index)
## table.gls.ts.atc<-apc.fit.table.2s(data.list.f,data.list.m,"gls.log.normal.rates","difference","APC","ATC",c(s.f/s.m,1),time.series=F_index)
## table.ols.ts.atc<-apc.fit.table.2s(data.list.f,data.list.m,"log.normal.rates","difference","APC","ATC",time.series=F_index)


###################################################
### code chunk number 19: ReproducingN2025.Rnw:325-333 (eval = FALSE)
###################################################
## # Combine analysis of variance results
## Table5<-cbind(table.ols.ts.apc[1:3,c(1,2,4,6,7)],
## rbind(c(NaN,NaN,NaN),
## table.ols.ts.atc[1:2,c(4,6,7)]),
## table.ols.ts.apc[1:3,10],
## table.gls.ts.apc[1:3,10])
## colnames(Table5)[c(1,9:10)]<-c("-2logL_OLS","sigma_OLS","sigma_GLS")
## Table5