% File apc/vignettes/CheckPrediction.Rnw % Part of the apc package, http://www.R-project.org % Copyright 2025 Bent Nielsen % Distributed under GPL 3 %\VignetteIndexEntry{Reproducing N2025} \documentclass[a4paper,twoside,12pt]{article} \usepackage[english]{babel} \usepackage{booktabs,rotating,graphicx,amsmath,verbatim,fancyhdr,Sweave} \usepackage[colorlinks,linkcolor=red,urlcolor=blue]{hyperref} \newcommand{\R}{\textsf{\bf R}} \renewcommand{\topfraction}{0.95} \renewcommand{\bottomfraction}{0.95} \renewcommand{\textfraction}{0.1} \renewcommand{\floatpagefraction}{0.9} \DeclareGraphicsExtensions{.pdf,.jpg} \setcounter{secnumdepth}{1} \setcounter{tocdepth}{1} \oddsidemargin 1mm \evensidemargin 1mm \textwidth 160mm \textheight 230mm \topmargin -5mm \headheight 8mm \headsep 5mm \footskip 15mm \begin{document} \SweaveOpts{concordance=TRUE} \raggedleft \pagestyle{empty} \vspace*{0.1\textheight} \Huge {\bf Reproducing \\ Nielsen (2022)\\ using the \texttt{apc} package} \noindent\rule[-1ex]{\textwidth}{5pt}\\[2.5ex] \Large 3 July 2025 \vfill \normalsize \begin{tabular}{rl} Bent Nielsen & Department of Economics, University of Oxford \\ & \small \& Nuffield College \\ & \normalsize \texttt{bent.nielsen@nuffield.ox.ac.uk} \\ & \url{http://users.ox.ac.uk/~nuff0078} \end{tabular} \normalsize \newpage \raggedright \parindent 3ex \parskip 0ex \tableofcontents \cleardoublepage \setcounter{page}{1} \pagestyle{fancy} \renewcommand{\sectionmark}[1]{\markboth{\thesection #1}{\thesection \ #1}} \fancyhead[OL,ER]{\sl Reproducing Nielsen (2025).} %\fancyhead[ER]{\sl \rightmark} \fancyhead[EL,OR]{\bf \thepage} \fancyfoot{} \renewcommand{\headrulewidth}{0.1pt} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} The purpose of this vignette is to use the \texttt{apc} package version 3.0.0 to reproduce some the result in Nielsen (2022): \textit{Two-sample age-period-cohort models with an application to Swiss suicide rates.}. The \texttt{apc} package builds on the identification analysis and the forecast theory in Kuang, Nielsen and Nielsen (2008a,b), the development of deviance analysis for general data arrays in Nielsen (2014). The package is discussed in Nielsen (2015). Fannon and Nielsen (2019) provide a review. The data originates from Riebler, Held, Rue and Bopp (2012). The data consists of counts of Swiss suicides and population by sex, in 5 year age groups from 15 to 79, by period from 1950 to 2007. It also includes time series such as a family index, the F-index. The purpose is to analyze if the age-period-cohort structure is common for the sexes and if the period effect can be replaced by a time series. This is modelled using a dose-response normal model with a two-sample age-period-cohort structure. The data are available in the \texttt{apc} package. They can be called with the command <<>>= # 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 @ Here \texttt{data.Swiss.suicides()} is a function that returns a list includes two \texttt{apc.data.list}s for females and males as well as five \texttt{vector}s of time series. Each \texttt{apc.data.list} is a list in itself including matrices for dose and response in \texttt{APm} format. This is a mixed scale format with $G=5$ year age groups $15-19$, $20-24$, $\dots$, $75-79$ and $H=1$ year periods $1950$, $1951$, $\dots$, $2007$. To see the structure of the function use the code <<>>= # Show some of the data objects(data) (data.list.f)$response[1:5,1:5] @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Figure 1: Crude rates per 100,000} The first figure shows crude suicide rates by age and period for women and men. The following command show plots separately for women and for men. <>= # Various plots of data sums apc.plot.data.sums(data.list.f,"c",scale.rate=100000) @ The figure in the paper has the crude rates for women and men in the same plot. First, a code for making two plots in one. <>= # 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 ) } @ Then, we compute the crude rates and scale with 100000. <<>>= # 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) @ And finally generate the plots <>= # Plot of data sums by age plot.both(crude.age.f,crude.age.m,"Crude rate by age","topleft") @ <>= # Plot of data sums by period plot.both(crude.per.f,crude.per.m,"Crude rate by period","topright") @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Table 2: Analysis of variance, women} <<>>= # 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 @ <<>>= # 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),] @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Table 3: Analysis of variance, men} The analysis of variance for men. The first command gives the basic table. The few lines add the estimated residual variances. <<>>= # 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 @ <<>>= # 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),] @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Table 4: Analysis of variance, both} The analysis of variance for both samples. First, it is done with OLS, that is imposing a common variance for the two samples. Second, it is done with GLS, that allows for different variance for the two samples. This requires that one specifies the relative variance. For this, the variances are computed from two one-sample analyses. Then the tables are put together. <>= # 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 @ Normality tests are computed for the two samples. <>= # 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),] @ \begin{table}[tb] \centering \begin{tabular}{lrrrrrrr} \hline \multicolumn{3}{l}{Model} & \multicolumn{5}{l}{Test vs APC} \\ diff.\ & $-2\log L_{GLS}$ & df & F & df & $p_F$ & $\hat\sigma_{OLS}$ & $\hat\sigma_{GLS}$ \\ % model & $-2\log L_{OLS}$ & df & F vs.\ apc & df vs.\ apc & $p_F$ & $\hat\sigma_{OLS}$ & $\hat\sigma_{GLS}$ \\ \hline APC & -1129.99 & 1144 & & & & 0.191 & 0.159 \\ AP & -836.72 & 1256 & 2.19 & 112 & 0.0000 & 0.201 & 0.167 \\ AC & -1052.62 & 1196 & 1.16 & 52 & 0.2092 & 0.192 & 0.160 \\ PC & -1021.68 & 1155 & 7.74 & 11 & 0.0000 & 0.197 & 0.164 \\ \hline \multicolumn{8}{l}{OLS: $\chi^2_{normality}(2)=15.82$ $[p=0.0004]$ } \\ \multicolumn{8}{l}{GLS: $\chi^2_{normality}(2)=\hspace{6pt}2.13$ $[p=0.3440]$ } \\ \hline \end{tabular} \caption{Two-sample analysis of variance. Submodels refer to restrictions on cross-sample differenced predictor, while the common predictor is an unrestricted APC model. For GLS, data for women are scaled to have same dispersion as men. \\ The structure of the table is close to that of Table \ref{tab:deviance:women}. Column 7 shows the scale estimated by OLS, so that the scale is common across samples. Column 8 shows the scale estimated by GLS, allowing different scales for the samples.} \label{tab:deviance:both} \end{table} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Figure 2: Detrended macro effects} This figure is generated from the two one sample fits. First, we plot the macro effects. At first, the effects from both samples are plotted together. There is one plot for each of age, period and cohort. Next, the cross-sample difference is plotted. Again, there are separate plots for each of age, period and cohort. If the option which.plot is omitted or set to 0 all six plots are made in one go. Otherwise, individual plots can be picked by choosing a value in the range 1-6 for which.plot. <>= # 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) @ Here the macro effects are identified by setting first and last value to zero. This preserves the degrees of freedom from the underlying double differences and ensures that there are no cross-constraints between plots. This also supports the interpretation of the time effects as showing non-linear effects only. Thus deviations from zero and the shape are interpretable. In line with this, the error bands are plus/minus to standard errors around zero. An alternative identification scheme is to set the two last age values and the first two period and cohort values to zero. This shows the underlying double sums of double difference. This is achieved by setting type to macro.ssdd. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Figure 3: Detrended micro effects} Next, we plot micro effects. The micro effects only exist when the data have mixed frequency as here. With $G=5$ age groups $H=1$ period groups only period and cohort micro effects exists. <>= # 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) @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Figure 4: F index} We will seek to replace the period effect with the F index. Only the non-linear part of the period effect is identified. Therefore only the non-linear part of the F index is relevant. Thus, we will plot a detrended version of the F index. We start by detrending <<>>= # Detrending F_index x.F <- 1950:2007 F_index.detrend <- F_index-F_index[1] - (0:57)/57*(F_index[58]-F_index[1]) @ We now plot the original F index along with the detrended version. The detrended version is multiplied by -1 to match macro period effect better. <>= # 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)) @ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Table 5: Analysis of variance with external time series} We do an analysis of variance where the period effect is replaced by an external time series. This is done by first computing four analyses of variance. We estimate by OLS and by GLS. Note that F values are the same with the two estimation methods. <>= # 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) @ We then combine the results. <>= # 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 @ \begin{table}[tb] \centering \begin{tabular}{lrrrrrrrrrr} \hline \multicolumn{3}{l}{Model} & \multicolumn{3}{l}{Test vs APC} & \multicolumn{5}{l}{Test vs AC+F} \\ diff.\ & $-2\log L_{GLS}$ & df & F & df & $p_F$ & F & df & $p_F$ & $\hat\sigma_{OLS}$ & $\hat\sigma_{GLS}$ \\ \hline APC &-1129.99 &1144 & & & & & & &0.191 &0.159 \\ AC+F &-1058.80 &1195 &51 &1.084 &0.321 & & & &0.191 &0.159 \\ AC &-1052.62 &1196 &52 &1.158 &0.209 &1 &4.906 &0.027 &0.192 &0.160 \\ \hline \end{tabular} \caption{Two-sample analysis of variance. Submodels refer to restrictions on cross-sample differenced predictor, while the common predictor is an unrestricted APC model. The AC+F model has the cross-sample difference period effect replaced with F-index. The structure of the table is close to that of Table \ref{tab:deviance:both}. Columns 4-6 report F-tests against the unrestricted APC model. Columns 7-9 report F-tests against the AC+F model.} % Submodels % refer to restrictions on cross-sample differenced predictor, while % the common predictor is an unrestricted APC model. % For GLS, data for women are scaled to have same dispersion as men.} \label{tab:deviance:external} \end{table} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section*{References} \begin{description} \item Fannon, Z. and Nielsen, B. (2019) Age-period-cohort models. Oxford Research Encyclopedia of Economics and Finance. Oxford University Press. Download: Article: \url{https://www.tandfonline.com/doi/full/10.1080/01621459.2017.1366908}. Early version circulated as Nuffield Discussion Paper 2018-W04 \url{https://www.nuffield.ox.ac.uk/economics/Papers/2018/2018W04_age_period_cohort_models.pdf}. \item Kuang, D., Nielsen, B. and Nielsen, J.P. (2008a) Identification of the age-period-cohort model and the extended chain ladder model. \textit{Biometrika} 95, 979-986. \textit{Download}: Earlier version: \url{http://www.nuffield.ox.ac.uk/economics/papers/2007/w5/KuangNielsenNielsen07.pdf}. \item Kuang, D., Nielsen, B. and Nielsen, J.P. (2008b) Forecasting with the age-period-cohort model and the extended chain-ladder model. \textit{Biometrika} 95, 987-991. \textit{Download}: Earlier version: \url{http://www.nuffield.ox.ac.uk/economics/papers/2008/w9/KuangNielsenNielsen_Forecast.pdf}. \item Nielsen, B.\ (2014) Deviance analysis of age-period-cohort models. \textit{Download}: \url{http://www.nuffield.ox.ac.uk/economics/papers/2014/apc_deviance.pdf}. \item Nielsen, B.\ (2015) apc: An R package for age-period-cohort analysis. \textit{R Journal} 7, 52-64. \textit{Download}: \url{https://journal.r-project.org/archive/2015-2/nielsen.pdf}. \item Nielsen, B. (2022) Two-sample age-period-cohort models with an application to Swiss suicide rates. \textit{Download}: Nuffield Discussion Paper 2022-W03: \url{https://www.nuffield.ox.ac.uk/economics/Papers/2022/2022-W03apc_2sample.pdf}. \item Riebler, A. and Held, L. and Rue, H. and Bopp, M.\ (2012) Gender-specific differences and the impact of family integration on time trends in age-stratified {S}wiss suicide rates. \textit{Journal of the Royal Statistical Society, Series A}, 175, 479-490. \end{description} \end{document}