\SweaveOpts{results=verbatim,keep.source=TRUE,include=FALSE,eps=FALSE,prefix.string=./prsz-}
\documentclass[a4paper,twoside,12pt]{article}

% ----------------------------------------------------------------------
% General information for the title page and the page headings
\newcommand{\Title}{Sample size and precision in the Greenland
  complications study}
\newcommand{\Tit}{Greenland precision}
\newcommand{\Version}{1}
\newcommand{\Dates}{October 2014}
\newcommand{\Where}{SDC / MaEJ}
\newcommand{\Homepage}{\url{http://bendixcarstensen.com/SDC/Greenland} }
\newcommand{\Faculty}{\begin{tabular}{rl}
Bendix Carstensen
  & Steno Diabetes Center, Gentofte, Denmark\\
  & {\small \& Department of Biostatistics,
               University of Copenhagen} \\
  & \texttt{bxc@steno.dk}\\
  & \url{http://BendixCarstensen.com} \\[1em]
                      \end{tabular}}

\input{/home/bendix/util/tex/toparticle}
\renewcommand{\rwpre}{./prsz-}

%----------------------------------------------------------------------
% Here comes the substance part(s)
%
\section{The study}

The Greenland complications study is an extra follow-up of the B99 and the
IHIT cohorts, where persons with diabetes or impaired glucose
regulation (IGR) are followed up for microvascular complications, notably
nephropathy and retinopathy.

The number of persons in the IGR is about 1200 (not known till the 2014
survey is completed), with an estimated follow-up time of 7,940
person-years \emph{after} onset of DM/IGR.

\section{Outcomes}

We want to estimate:
\begin{itemize}
\item rates of microvascular complications
\item rate-ratios of microvascular complications between carriers of different
 genotypes of interest.
\end{itemize}

\subsection{Assumptions}

Based on previous studies and sound clinical judgment the rate of
newly onset nephropathy among persons with DM or IGR is taken to be 20
per 1000 PY and the corresponding rate of retinopathy to be 30 per
1000 PY.

\section{Precision and power}

\subsection{Rates}

A bold calculation yields that a complication with a rate of 20/1000
during 7,940 PY will produce $20 \times 7.94 = 159$ events, and one
with rates 30/1000 will produce 238 events:
<<>>=
c(20,30) * 7.94
@ %
The usual formula for the relative precision of a rate estimate based on
these numbers is:
\[
  \erf = \exp( 1.96/\sqrt{D} )
\]
which in the cases here will produce rate estimates with a relative
precision in the vicinity of $\pm$ 15\%:
<<>>=
exp( 1.96/sqrt(c(159,238)) )
@ %

\subsection{Rate-ratios}

We are looking at genotypes with population prevalences between 5 and
10\% separately and 44\% in total, but we do not have any specific
estimates of the likely risk elevation associated each of them.

Hence we will compute the precision of the RR estimates using a
baseline rates of 20 and 30 per 1000 PY, for exposure prevalences 5 --
50 \%, and exposure effects 1.05, \ldots, 1.50. This will slightly
overestimate the precision if the true overall rates is 20 resp. 30,
because a fraction of the population having an elevated risk will of
course result in an overall rate that is higher than the rate in the
unexposed group. For the sake of completeness we also make the
computations for person-years between 8,000 and 13,000.

Thus, we set up an array to collect the median RR (for control), the
precision (in terms of the error factor, the relative uncertainty of
the RR estimates) and the power for combinations of these factors:
<<>>=
library( Epi )
nsim <- 10000
prpw <- NArray( list( RR = seq(1.05,1.5,0.05),
                      pr = seq(5,50,5),
                      PY = 8:13,
                     inc = c(20,30),
                    what = c("RR","Erf","Pwr") ) )
str( prpw )
@ %
With this array to collect the results we can simulate, by approximating
the no. events by a Poisson variate with mean equal to the cumulative rate:
<<>>=
for( ir in dimnames(prpw)[["RR"]] )
for( ip in dimnames(prpw)[["pr"]] )
for( iy in dimnames(prpw)[["PY"]] )
for( ii in dimnames(prpw)[["inc"]] )
{
nr <- as.numeric(ir)
np <- as.numeric(ip)/100
ny <- as.numeric(iy)
ni <- as.numeric(ii)
# Cumulative incidence in exposed and non-exposed:
cuminc <- ny * c(np,1-np) * ni * c(nr,1)
# Events in exposed:
nx <- rpois( nsim, cuminc[1] )
# Events in unexposed:
nr <- rpois( nsim, cuminc[2] )
# Rate Ratio
rr <-  (nx/np) / (nr/(1-np))
# Error factor
erf <- exp( 1.96*sqrt(1/nx+1/nr) )
# Significance
sgn <- abs(log(rr)) > log(erf)
# Store results
prpw[ir,ip,iy,ii,"RR" ] <- median( rr )
prpw[ir,ip,iy,ii,"Erf"] <- median( erf )
prpw[ir,ip,iy,ii,"Pwr"] <- mean( sgn )
}
@ %
For the sake of verification, we print the median RRs from the
simulations:
<<>>=
round( ftable( prpw[,,"12",,"RR" ],
               col.vars=1, row.vars=c(3,2) ), 2 )
@ %
Further more we see that the precision of the estimates (the
error-factor) decreases only little with increasing RR, but much more
so by increasing prevalence of the risk factor:
<<>>=
round( ftable( prpw[,,c("8","12"),,"Erf"],
               col.vars=1, row.vars=c(3,4,2) ), 2 )
@ %
We can then print the anticipated power to detect a significant RR of
a given size, in percent:
<<>>=
round( ftable( prpw[,,c("8","12"),,"Pwr"]*100,
               col.vars=1, row.vars=c(3,4,2) ), 1 )
@ %
And show where the power exceeds 80\%:
<<>>=
print( ftable( (prpw[,,c("8","12"),,"Pwr"]>0.8)*1,
               col.vars=1, row.vars=c(3,4,2) ),
       zero.print="." )
@ %

Broadly speaking, with 80\% power we can only see effects of RR$>1.4$
for exposures with prevalence over 30\%, but the good news is that
with 12,000 PY we have a relative precision of 1.3 and with 8,000 PY a
relative precision in the vicinity of 1.5.

% \bibliographystyle{plain}
% \bibliography{%
% c:/bendix/artikler/bibtex/BxC,%
% c:/bendix/artikler/bibtex/diabetes%
% }
% \addcontentsline{toc}{section}{References}

\end{document}