\SweaveOpts{results=verbatim,keep.source=TRUE,include=FALSE,eps=FALSE,prefix.string=./graph/CIMT} \documentclass[a4paper,twoside,12pt]{report} % ---------------------------------------------------------------------- % General information for the title page and the page headings \newcommand{\Title}{CIMT study\\ Ultrasound and metabolic data} \newcommand{\Tit}{CIMT} \newcommand{\Version}{1 (preliminary)} \newcommand{\Dates}{September 2013} \newcommand{\Where}{SDC / LLCh} \newcommand{\Homepage}{\url{http://BendixCarstensen.com/SDC/CIMT} } \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{c:/util/tex/topreport} \renewcommand{\rwpre}{./graph/CIMT} %---------------------------------------------------------------------- % Here comes the substance part(s) % \chapter{Construction of data} First we read the ultrasound data for \emph{all} participants. Note that some of the numerical variables are with ``.'' as decimal separator, even if we have a file from a Danish Locale: <<>>= imt <- read.csv2("./data/imt.csv") str( imt ) names( imt )[c(5,9,12)] for( i in c(5,9,12) ) imt[,i] <- as.numeric(as.character(imt[,i])) imt$datescanned <- as.Date( substr(imt$datescanned,1,10) ) str(imt) addmargins( with( imt, table(visit,cca_side) ) ) with( imt, table( table( subjid ) ) ) @ % The we read the randomization codes for those actually randomized: <<>>= rnd <- read.csv2("./data/PTNR_MetforminCode_InsulinTypeCode.csv" ) str( rnd ) @ % By merging with \texttt{all.y=TRUE}, we ensure that only those from the randomization set are included: <<>>= dim( imt ) rimt <- merge( imt, rnd, by.x="subjid", by.y="PTNR", all.y=TRUE ) dim( rimt ) str( rimt ) with( rimt, table(table(subjid)) ) @ % Fishily enough there is one person with 4 observations in \texttt{imt} that is not randomized: <<>>= tt <- with( imt, table( subjid ) ) has4i <- names( tt[tt==4] ) isrnd <- rnd$PTNR setdiff( has4i, isrnd ) @ % <>= with( rimt, pairs( cbind( avg= as.numeric(as.character(mean_fimtavg)), min=mean_fimtmin, max=mean_fimtmax, gns=(mean_fimtmin+mean_fimtmax)/2 ), gap=0, pch=16, cex=0.5 ) ) @ % $\ \insfig{pairs}{0.9}{Pairwise comparison of far wall intima media thickness as the average, min, max and the mean of min and max.} We now take the average for left and right-sided measurements, but also merge in as separate variables the left and right side measurements: <<>>= what <- 5:22 names( rimt )[what] mimt <- aggregate( rimt[,what], list( subjid=rimt$subjid, visit=rimt$visit), FUN=mean, na.rm=TRUE ) wh <- c(1,2,5:12,18:20) names( imt )[wh] Limt <- subset( rimt, cca_side=="L" )[,wh] Rimt <- subset( rimt, cca_side=="R" )[,wh] names( Limt ) <- gsub( "mean_", "", names( Limt ) ) names( Rimt ) <- gsub( "mean_", "", names( Rimt ) ) names( Limt )[-(1:2)] <- paste( "left_" , names( Limt )[-(1:2)], sep="" ) names( Rimt )[-(1:2)] <- paste( "right_", names( Rimt )[-(1:2)], sep="" ) mimt <- merge( merge( mimt, Limt, all.x=TRUE ), Rimt, all.x=TRUE ) head( mimt ) subset( rimt, subjid==10001 ) subset( mimt, subjid==10001 ) @ % Finally we attach the stratification variables: <<>>= st <- read.csv2( "./data/strata.csv" )[,-2] names( st )[2:4] <- c("o.65","pre.ins","SDC") mimt <- merge( mimt, st, by.x="subjid", by.y="ptnr" ) @ % The resulting data set has one observation per visit among the randomized patients: <<>>= with( mimt, ftable(addmargins(table(Met=MetforminCode, Ins=InsulinTypeCode, visit ))) ) @ % Finally for the benefit of the client, we write it both an \texttt{Rda} format and as a SAS data-file with corresponding script: <<>>= save( mimt, file="./data/mimt.Rda" ) write.foreign( mimt, datafile = "./sas/mimt.dat", codefile = "./sas/mimt.sas", package = "SAS" ) @ % $ \chapter{Primary outcome (CIMT)} \section{Analysis as repeated measures} We reload the data and the two necessary packages from \R, and convert the subject indicator to a factor as is needed for use in \texttt{lmer}: <<>>= library( lme4 ) library( Epi ) load( file="./data/mimt.Rda" ) mimt <- transform( mimt, subjid = factor(subjid), met = factor(MetforminCode,labels=c("Plc","Met")), ins = factor(InsulinTypeCode,labels=c("NR+Lv","Lv","NM30")) ) with( mimt, table(MetforminCode, met )) with( mimt, table(InsulinTypeCode,ins)) @ % We analyze data with a random effects model, for the CIMT-mean $y_{it}$ on individual $i$ at time $t=1,7$ (\texttt{mean\_fimtavg}), randomized to treatment $m=M,P$ using subject as random, and with a separate metformin by time (\texttt{visit}) interaction: \begin{align*} y_{it} & = \mu + \beta_t + \gamma_{mt} \\ & \quad + \alpha_1\mathtt{o.65} + \alpha_2\mathtt{per.ins} + \alpha_1\mathtt{SDC}\\ & \quad + a_i + e_{it},\\ & a_i \sim \mathcal{N}(0,\tau^2), \quad e_{it} \sim \mathcal{N}(0,\sigma^2), \quad \text{all independent} \end{align*} The model states that persons have an average baseline level of CIMT depending on randomization and stratification group; the mean at baseline and follow-up are: \begin{align*} \text{baseline} & \mu+\beta_1+\gamma_{m1}+\alpha_1\mathtt{o.65}+\alpha_2\mathtt{per.ins}+\alpha_3\mathtt{SDC} \\ \text{follow-up} & \mu+\beta_7+\gamma_{m7}+\alpha_1\mathtt{o.65}+\alpha_2\mathtt{per.ins}+\alpha_3\mathtt{SDC} \\ \text{change} & \beta_7-\beta_1 + (\gamma_{m7}-\gamma_{m1}) \end{align*} So in the changes in CIMT are: \begin{align*} \text{metformin} & \beta_7-\beta_1 + (\gamma_{\text{met}7}-\gamma_{\text{met}1}) \\ \text{placebo} & \beta_7-\beta_1 + (\gamma_{\text{plc}7}-\gamma_{\text{plc}1}) \end{align*} The model as stated here is grossly overparametrized, it can be identified if all parameters relation to baseline or placebo were set to 0, the changes in CIMT over the follow-up are \begin{align*} \text{metformin} & \beta_7 + \gamma_{\text{met}7} \\ \text{placebo} & \beta_7 \\ \text{difference} & \gamma_{\text{met}7} \end{align*} These three parameters are those of interest which should be extracted. The model is a random effects model that is very close to using the baseline ($y_{i1}$) as covariate in an analysis of the follow-up ($y_{i7}$) as outcome: <<>>= m0 <- lmer( mean_fimtavg ~ -1 + visit:met + o.65 + pre.ins + SDC +(1|subjid), data = mimt ) summary( m0 ) round( ee <- ci.lin( m0, subset="visit" ), 4 ) CM <- rbind(diag(4),c(0,0,-1,1), c(-1,1,0,0), c(1,-1,-1,1)) row.names(CM) <- c("Plc 1a","Plc 7a", "Met 1a","Met 7a", "Met 7a-1a", "Plc 7a-1a", "dMet - dPlc" ) colnames(CM) <- row.names(ee) CM round( e0 <- ci.lin( m0, subset="visit", ctr.mat=CM ), 4 ) elab <- c("Metformin 18m - baseline", "Placebo 18m - baseline", "Metformin vs. Placebo change") rownames( e0 )[5:7] <- elab round( e0, 4 ) @ % <<>>= eclr <- c("green","yellow","orange") win.metafile( "forest1.emf", width=10, height=5, pointsize=24 ) par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6, bg="black","transparent", col.axis="white", col.lab="white" ) plotEst( e0[5:7,c(1,5,6)], lwd=7,cex=1.5, xlab="", col.txt="transparent", col=eclr, restore.par=FALSE ) for( i in 1:3 ) axis( side=2, at=4-i, labels=elab[i], col="transparent", col.axis= eclr[i], las=1 ) abline(v=0,col="white") axis( side=1, col="white" ) mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6, col="white" ) dev.off() @ <>= par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 ) plotEst( e0[5:7,c(1,5,6)], lwd=3, cex=1.5, vref=0, xlab="Carotid intima-media thickness (mm)", txt=elab ) @ % \insfig{forest1}{0.9}{Estimated contrasts from the random effects model.} \section{Analysis using baseline as covariate} Another way of looking at the analysis is to use the measurement at the second occasion using the baseline measurement as covariate. In section \ref{sex:relation} is a brief account of the relationship between the two approaches. \subsection{Model} The main difference is that the absolute change in the response from baseline to follow-up in each of the groups does not appear as a parameter in this approach. The model is (for $F$---follow-up and $B$---baseline): \begin{align*} y_{iF} & = \mu + \theta y_{iB} + \delta_t \\ & \quad + \alpha_1\mathtt{o.65} + \alpha_2\mathtt{per.ins} + \alpha_1\mathtt{SDC}\\ & \quad + e_i,\\ & e_i \sim \mathcal{N}(0,\sigma^2) \end{align*} % Thus in the placebo group the change in $y$ is $y_{iF}-y_{iB}$, which % in this model equals $(\theta-1)y_{iB}$, so not an absolute but a % \emph{relative} change by a factor $(\theta-1)$, whereas the change in % the treatment group is $(\theta-1)y_{iB}+\delta_t$. So the model % assumes that the \emph{relative} change from baseline to follow-up is % constant across individuals, and that the treatment effect is % \emph{additive} on top of this. To the extent that other variables are % included in the model these model the additive effect on the change too. % Hence there is no such thing as the absolute change in each of the % groups in this model. In order to use this approach we restructure the data set so that we have only one row per person with follow-up and baseline as separate variables: <<>>= names( mimt ) wimt <- reshape( mimt[,c("mean_fimtavg","visit","subjid", "MetforminCode","o.65","pre.ins","SDC")], direction = "wide", v.names = "mean_fimtavg", timevar = "visit", idvar = "subjid" ) subset( mimt[,c("mean_fimtavg","visit","subjid", "MetforminCode","o.65","pre.ins","SDC")], subjid=="10002" ) subset( wimt, subjid=="10002" ) @ % We now have a dataset with the relevant variables, where we can estimate the difference in changes between the two groups: <<>>= mf <- lm( mean_fimtavg.7a ~ mean_fimtavg.1a + MetforminCode + o.65 + pre.ins + SDC, data = wimt ) summary( mf ) mm <- lm( mean_fimtavg.7a ~ mean_fimtavg.1a + MetforminCode, data = wimt ) summary( mm ) round( cf <- ci.lin( mf, subset=c("fimt","Met") ), 4 ) round( cm <- ci.lin( mm, subset=c("fimt","Met") ), 4 ) @ % These estimates can now be compared with those from the random effects model: <<>>= round( e0, 4 ) @ % The relevant quantities to compare are the additive effects of treatment: <<>>= round( rbind( "Simple"=cf[2,], "Contrl"=cm[2,], "RanEff"=e0[7,] ), 4 ) @ % Clearly, the conclusion for the models are not substantially different; CIMT difference is some 0.01 mm larger between the groups at follow-up in favor of the Placebo group. To understand the estimates pertaining to the baseline variable, we can compare them to the ratio of the main effects from the linear mixed model. % Note however that the ratio of follow-up to baseline in % the two treatment groups is different, whereas it in the simplfied % model is assumed constant. \section{Using multiple imputation} We see that there are some missing values in the follow-up values. Therefore we want to make a multiple imputation using a rich model to predict the outcome for additional analyses. Hence we are using the multiple imputation <>= options( width=90 ) summary( wimt ) @ We start by reading the data set with all the variables to be used for imputation, except the outcomes variables CIMT: <<>>= library( mice ) library( foreign ) llch <- read.ssd( "./data", "mi" ) names( llch ) <- tolower( gsub("_",".",names(llch)) ) mimt <- merge( wimt, llch ) names( mimt ) <- gsub( "mean_", "", names(mimt) ) mimt <- transform( mimt, sex = factor( sex , labels=c("F","M") ), # code 0/1 ! ami = factor( ami , labels=c("No","Yes") ), heartsur = factor( heartsur, labels=c("No","Yes") ), apopl = factor( apopl , labels=c("No","Yes") ), tci = factor( tci , labels=c("No","Yes") ), sdc = factor( sdc , labels=c("no","yes") ), o.65 = factor( o.65 , labels=c("no","yes") ), pre.ins = factor( pre.ins , labels=c("no","yes") ), metformi = factor( metformi, labels=c("Plc","Met") ) ) mimt <- mimt[,c(1,42,43,41,3,4,6:40,44:47)] cbind( 1:ncol(mimt), t( apply( mimt, 2, function(x) c(sum(is.na(x)), length(unique(x))) ) ) ) ftable( addmargins( xtabs( cbind( "Baseline"=!is.na(fimtavg.1a), "FU18mth"=!is.na(fimtavg.7a), "Diff."=!is.na(fimtavg.1a)- !is.na(fimtavg.7a), "MissingFU"= is.na(fimtavg.7a) ) ~ sex + metformi, data=mimt ), margin=1:2 ) ) @ % So we see that there is about 10\% missing follow-up measurements (41 out of 412). The pattern of missing values is: <<>>= zz <- md.pattern( mimt ) t( zz ) @ <>= pairs( mimt[,c(7,8,15,16:24)], gap=0, pch=16, cex=0.2 ) @ % <>= pairs( mimt[,25:36], gap=0, pch=16, cex=0.2 ) @ % <>= pairs( mimt[,37:45], gap=0, pch=16, cex=0.2 ) @ % \insfig{pairs1}{1.0}{Pairwise marginal distributions of the variables in the dataset \textrm{\tt mimt}.} \insfig{pairs2}{1.0}{Pairwise marginal distributions of the variables in the dataset \textrm{\tt mimt}.} \insfig{pairs3}{0.75}{Pairwise marginal distributions of the variables in the dataset \textrm{\tt mimt}.} We see from figure \ref{fig:pairs1} that non-severe hypoglycemia, \texttt{nonsever} is quite skew so we will include it as log-transformed in the imputation. In figure \ref{fig:pairs2} it is evident that the weight measurements very closely correlated, so we include only the baseline and the visit7 values. Finally from figure \ref{fig:pairs3} we see that the albumin measurements are both correlated and very skewed, so we include only the log-transform of the average (\texttt{avgnatua}). With this in the model we do not include microalbuminuria: <>= mimimt <- transform( mimt, l.alb = log( avgnatua ), l.nsev = log( nonsever ) )[,-c(14,16,26:30,42:45)][c(1:15,36,16:35)] names( mimimt ) save( mimimt, file="./data/mimimt.Rda" ) @ % <>= pairs( mimimt[,7:21], gap=0, pch="." ) @ % <>= pairs( mimimt[,22:36], gap=0, pch="." ) @ % With this in place we can now make a multiple imputation of the missing values, including those on the outcome variable: <>= set.seed( 876528358 ) system.time( imp <- mice( mimimt[,2:36], m=50 ) ) save( imp, file="./data/imp" ) @ % Once we have done the imputation we can make the imputation analysis: <>= class( imp ) i.mf <- with( imp, lm( fimtavg.7a ~ fimtavg.1a + metformi + o.65 + pre.ins + sdc ) ) round( smf <- summary( pool(i.mf) ), 4 ) i.mm <- with( imp, lm( fimtavg.7a ~ fimtavg.1a + metformi ) ) round( smm <- summary( pool( i.mm ) ), 4 ) @ With the estimate for the metformin effect in \texttt{smm} (which is the relevant one) we can now re-do the forest plot using the ``old'' estimates of the corrected changes estimated from the random-effects model: <>= par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6 ) plotEst( rbind(e0[5:6,c(1,5,6)],smm[3,c("est","lo 95","hi 95")]), txt=rownames(e0)[5:7], lwd=7,cex=1.5, restore.par=FALSE ) abline( v=0 ) axis( side=1 ) mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6 ) @ \insfig{forest2}{0.9}{Estimated contrasts from the random effects model. The difference is the estimate from the model with baseline as covariate, using multiply imputed data.} <<>>= eclr <- c("green","yellow","orange") win.metafile( "forest2.emf", width=10, height=6, pointsize=20 ) par( mar=c(3,1,1,1), mgp=c(3,1,0)/1.6, bg="transparent", col.axis="white", col.lab="white" ) plotEst( rbind(e0[5:6,c(1,5,6)],smm[3,c("est","lo 95","hi 95")]), ylim=c(0,3), txt=rownames(e0)[5:7], lwd=7,cex=1.5, xlab="", col.txt="transparent", col=eclr, restore.par=FALSE ) for( i in 1:3 ) axis( side=2, at=4-i, labels=elab[i], col="transparent", col.axis= eclr[i], las=1 ) abline(v=0,col="white") axis( side=1, col="white" ) text( c(e0[5:6,1],smm[3,"est"]), 0.7+2:0, paste( "P =", formatC( c(e0[5:6,"P"],smm[3,"Pr(>|t|)"]), format="f", digits=3 ) ), col="white" ) text( -0.03, 0.2, "Improvement" , adj=0, col="white", cex=0.8 ) text( 0.03, 0.2, "Deterioration", adj=1, col="white", cex=0.8 ) mtext( "Carotid intima-media thickness (mm)", side=1, line=3/1.6, col="white" ) dev.off() @ % \section{Relationship between between ``repeated measures'' and % ``baseline as covariate'' approaches} % \label{sec:relation} % As noted above the two approaches are substantially different, in that % one models the % If we take the random effects model as strating point, \chapter{Glucose-related variables} Here is the documentation of the reading of the data: \insout{./sas/gethbmv.yt} \section{Hb$_\text{A1c}$} <<>>= options( width=110 ) library( Epi ) library( foreign ) library( lme4 ) hba <- read.xport( "./data/hba.xpt" ) names( hba ) <- tolower( names( hba ) ) names( hba ) <- gsub( "_", ".", names( hba ) ) hba$grp <- factor( hba$metformi, labels=c("plc","met") ) names( hba )[2] <- "inspw" str( hba ) summary( hba ) @ % How many persons with how many measurements <<>>= ( tt <- with( hba, table(table(subjid)) ) ) round( 100*tt/sum(tt), 1 ) @ \subsection{The complete cases} We select for a start those with complete follow-up <<>>= vv <- with( hba, table(subjid) ) cplv <- names( vv[vv==7] ) hbc <- subset( hba, subjid %in% cplv ) with( hba, table(table(subjid)) ) with( hbc, table(table(subjid)) ) names( hbc ) @ We set up a random effects model using \texttt{lmer} <>= mc <- lmer( hba1c ~ grp + grp:factor(visit) - 1 + factor(sdc) + factor(above.si) + factor(prior.in) + (1|subjid), data=hbc ) ma <- update( mc, data=hba ) ( eM <- cbind(1,rbind(0,diag(6))) ) ci.lin( mc ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(-eM,eM) ) round( c.dif, 3 ) ci.lin( ma ) a.plc <- ci.lin( ma, ctr.mat=cbind(1,rbind(0,diag(6))), subset="plc" )[,-(2:4)] a.met <- ci.lin( ma, ctr.mat=cbind(1,rbind(0,diag(6))), subset="met" )[,-(2:4)] a.dif <- ci.lin( ma, subset=c("met","plc"), ctr.mat=cbind(-eM,eM) ) round( a.dif, 3 ) P.all <- Wald( ma, H0=rep(0,7), subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) P.last<- Wald( ma, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[7,,drop=F] ) win.metafile( "hba1c.emf", width=10, height=6, pointsize=20 ) eclr <- c("yellow","green","orange") par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:6*3, cbind( a.plc, a.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,19.5), ylim=c(7,9.5), yaxs="i", yaxt="n", ylab=expression( "Hb"*A[1][c]*" (%) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) # abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:6*3, cbind( a.plc, a.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=7:9, col="white" ) text( c(20,20), c(9.1,9.3), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=1, cex=1 ) text( 20, 8.8, substitute( "Equal Hb"*A[1][c]*" 0-18 mth: P ="*pval, list(pval = formatC( P.all[3], format="f", digits=3 )) ), col="white", adj=1 ) text( 20, 8.6,substitute( "Equal Hb"*A[1][c]*" at end: P ="*pval, list(pval = formatC( P.last[3], format="f", digits=3 )) ), col="white", adj=1 ) # text( 19.2, 8.55, "HbA1c change\n from baseline", col="white", adj=c(1,1) ) text( c(20,20), c(a.plc[7,1],a.met[7,1]), paste( formatC( c(a.plc[7,1],a.met[7,1])-c(a.plc[1,1],a.met[1,1]), format="f", digits=1 ), "%", sep="" ), col=eclr[1:2], adj=1 ) dev.off() @ % \section{Insulin dose } \subsection{The complete cases} We select for a start those with complete follow-up <<>>= vv <- with( hba, table(subjid) ) cplv <- names( vv[vv==7] ) hbc <- subset( hba, subjid %in% cplv ) with( hba, table(table(subjid)) ) with( hbc, table(table(subjid)) ) names( hbc ) with( hba, tapply( inspw, visit, function(x) sum( !is.na(x) ) ) ) with( hba, tapply( inspw, visit, function(x) sum( is.na(x) ) ) ) @ We set up a random effects model using \texttt{lmer} <>= mc <- lmer( inspw ~ grp + grp:factor(visit) - 1 + factor(sdc) + factor(above.si) + factor(prior.in) + (1|subjid), data = subset( hba, visit!=7) ) round( ci.lin( mc ), 3 ) ( eM <- cbind(1,rbind(0,diag(5))) ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) round( c.dif, 3 ) P.all <- Wald( mc, H0=rep(0,6), subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) P.last<- Wald( mc, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[6,,drop=F] ) @ % For annotation we compute the empirical mean insulin doses in each of the randomization groups, at visit 6: <<>>= ( IL <- with( hba, tapply( insdose2, list( metformi, visit), mean ) )[,6] ) @ % <<>>= eclr <- c("yellow","green","orange") win.metafile( "inspw.emf", width=10, height=6, pointsize=20 ) par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:5*3, cbind( c.plc, c.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,18), ylim=c(0,1.3), yaxs="i", yaxt="n", ylab=expression( "Insulin (IU/kg/day) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:5*3, cbind( c.plc, c.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=seq(0,1.2,0.2), col="white" ) text( c(18,18), c(0.35,0.45), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=1, cex=1 ) text( 18, 0.2, paste( "Equal insulin dose 0-15 mth: P =", formatC( P.all[3], format="f", digits=3 ) ), col="white", adj=1 ) text( 18, 0.1, paste( "Equal insulin dose at end: P =", formatC( P.last[3], format="f", digits=3 ) ), col="white", adj=1 ) text( c(18,18), c(c.plc[6,1],c.met[6,1]), paste( round(IL), "IU/day"), col=eclr[1:2], adj=1 ) dev.off() @ % \section{Weight} We use the same dataset as before: <<>>= names( hba ) with( hba, tapply( weight, visit, function(x) sum( !is.na(x) ) ) ) with( hba, tapply( weight, visit, function(x) sum( is.na(x) ) ) ) @ \subsection{Absolute weight} We set up a random effects model using \texttt{lmer} <>= mc <- lmer( weight ~ grp + grp:factor(visit) - 1 + factor(sdc) + factor(above.si) + factor(prior.in) + (1|subjid), data = hba ) round( ci.lin( mc ), 3 ) ( eM <- cbind(1,rbind(0,diag(6))) ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) round( c.dif, 3 ) P.all <- Wald( mc, H0=rep(0,7), subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) P.last<- Wald( mc, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[7,,drop=F] ) @ % <<>>= win.metafile( "weight.emf", width=10, height=6, pointsize=24 ) eclr <- c("yellow","green","orange") par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:6*3, cbind( c.plc, c.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,18), ylim=c(90,110), yaxs="i", yaxt="n", ylab=expression( "Weight (kg) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) # abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:6*3, cbind( c.plc, c.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=seq(90,110,5), col="white" ) text( c(17,17), c(0.25,0.35), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=1, cex=1 ) text( 17, 0.15, paste( "Equal insulin dose 0-15 mth: P =", formatC( P.all[3], format="f", digits=3 ) ), col="white", adj=1 ) text( 17, 0.1, paste( "Equal insulin dose at end: P =", formatC( P.last[3], format="f", digits=3 ) ), col="white", adj=1 ) dev.off() @ % \subsection{Weight change} The fitted model also contains the estimates of the weight change, it is merely a question of changing the contrast matrix, by replacing the 1-row that fishes out the baseline at time 0 by a row of 0s that fixes the effect at time 0 to 0, and hence the estimates to changes: <<>>= round( ci.lin( mc ), 3 ) ( eM <- cbind(0,rbind(0,diag(6))) ) c.plc <- ci.lin( mc, ctr.mat=eM, subset="plc" )[,-(2:4)] c.met <- ci.lin( mc, ctr.mat=eM, subset="met" )[,-(2:4)] c.dif <- ci.lin( mc, subset=c("met","plc"), ctr.mat=cbind(eM,-eM) ) round( c.dif, 3 ) P.all <- Wald( mc, H0=rep(0,6), subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[-1,] ) P.last<- Wald( mc, H0=0, subset=c("met","plc"), ctr.mat=cbind(eM,-eM)[7,,drop=F] ) @ <<>>= win.metafile( "wchg.emf", width=10, height=6, pointsize=20 ) eclr <- c("yellow","green","orange") par( mar=c(3,3,1,1), mgp=c(3,1,0)/1.6, las=1, bg=#"black", "transparent", col.axis="white", col.lab="white", bty="n" ) matplot( 0:6*3, cbind( c.plc, c.met ), xlab="Months since trial entry", xaxt="n", xlim=c(0,18), ylim=c(0,6.2), yaxt="n", ylab=expression( "Weight change (kg) [Mean (95% CI)]" ), type="n", lwd=c(4,1,1), lty=1, col=rep(eclr[1:2],each=3) ) # abline( h=seq(7,9.5,0.5), col=gray(0.5) ) matlines( 0:6*3, cbind( c.plc, c.met ), type="l", lwd=c(5,1,1), lty=1, col=rep(eclr[1:2],each=3) ) axis( side=1, at=0:6*3, col="white" ) axis( side=2, at=0:6, col="white" ) text( c(0.5,0.5), c(5.7,6.2), c("Placebo","Metformin"), col=eclr[1:2], font=2, adj=0, cex=1 ) text( 0.5, 5.2, paste( "Equal weight change 0-18 mth: P =", formatC( P.all[3], format="f", digits=3 ) ), col="white", adj=0 ) text( 0.5, 4.7, paste( "Equal weight change end: P =", formatC( P.last[3], format="f", digits=3 ) ), col="white", adj=0 ) dev.off() @ % \end{document} \bibliographystyle{plain} \bibliography{% c:/bendix/artikler/bibtex/BxC,% c:/bendix/artikler/bibtex/diabetes% } \addcontentsline{toc}{section}{References} \end{document}