diff --git a/DESCRIPTION b/DESCRIPTION
index c7b358d..6e8a8a6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -2,16 +2,16 @@ Package: TP.idm
 Type: Package
 Title: Estimation of Transition Probabilities for the Illness-Death
         Model
-Version: 1.2
-Date: 2016-11-28
+Version: 1.3
+Date: 2018-02-14
 Author: Vanesa Balboa-Barreiro, Jacobo de Una-Alvarez <jacobo@uvigo.es>
          and Luis Meira-Machado <lmachado@math.uminho.pt>
 Maintainer: Vanesa Balboa-Barreiro <vanesa.balboa@uvigo.es>
 Depends: R (>= 3.1.1), graphics, stats
 Repository: CRAN
-Description: Estimation of transition probabilities for the illness-death model. Both the Aalen-Johansen estimator for a Markov model and a novel non-Markovian estimator by de Una-Alvarez and Meira-Machado (2015) <DOI:10.1111/biom.12288> are included.
+Description: Estimation of transition probabilities for the illness-death model. Both the Aalen-Johansen estimator for a Markov model and a novel non-Markovian estimator by de Una-Alvarez and Meira-Machado (2015) <doi:10.1111/biom.12288>, see also Balboa and de Una-Alvarez (2018) <doi:10.18637/jss.v083.i10>, are included.
 License: GPL-2
 Encoding: latin1
 NeedsCompilation: yes
-Packaged: 2016-11-30 20:54:53 UTC; hp
-Date/Publication: 2016-12-01 08:28:20
+Packaged: 2018-02-19 20:25:16 UTC; hp
+Date/Publication: 2018-02-20 09:15:24 UTC
diff --git a/MD5 b/MD5
index dd9dc41..d18e61c 100644
--- a/MD5
+++ b/MD5
@@ -1,26 +1,29 @@
-84eebc8de2f233abeb8ca93a65dce83c *DESCRIPTION
-a983a957fbeca705bf96d3311974ecb9 *NAMESPACE
+54d2da93fe1e6558728f6eb137c97c03 *DESCRIPTION
+a6fb9a499ddf54c27a9cbaac5d0c92e9 *NAMESPACE
 b871681d0ffbb00da9ff7056f7577210 *R/KMW.R
-eb2736fb5cae22d84ff31ae67ad239dc *R/TPidm.R
-027079d0156146edd4940c4bf2d7e04d *R/ci.AJ.R
+0b03bd0ea45acd34e5842006ec4a9ca9 *R/TP.idm-internal.R
+0baa720622236970c34bbad967caefe9 *R/TPidm.R
+6788149f6617f0c6d7de2f0f8ba20ae7 *R/ci.AJ.R
 797154ed05e60313b5f389561e657db3 *R/ci.NM.R
-5e6a0851f4eba40677dc267e3169b72d *R/fun.AJ.R
+5d533496bc2f7b5b57fca9332283056f *R/fun.AJ.R
 8e9c27af6623eca95c2a00fb40b757f8 *R/fun.NM.R
-e039db4d105c6c0fc6ca228a43770d37 *R/plot.TPidm.R
-9885193d94b1989b826f7a9d39b2e239 *R/prep.data.AJ.R
-61515d3759c287cb1943b53dbb9b2c6d *R/prep.data.event.AJ.R
-db36af6c5cb01c07caf87d5a395560dd *R/print.TPidm.R
+9acb6fd878f73487f4b448aed6838857 *R/plot.TPidm.R
+336d078c9a471b31abe95800f385125c *R/prep.data.AJ.R
+ef5cd6f24ee725449c51afc9a864159e *R/prep.data.event.AJ.R
+69b9f171f646db2341fc5dfa69fb38e6 *R/print.TPidm.R
 d0a5320eeab15164381fa662c5674ab7 *R/summary.TPidm.R
-4d90cea95426cbfba2281ddc67e513b6 *R/test.nm.R
-b994c0a5a572e7420f8e1b6872b7ced2 *R/var.AJ.R
-51b24c0f0fc828bcf9bf29b1d056ddff *R/var.NM.R
-a3aa8f2a86ed57745e8b924837bead6b *build/partial.rdb
-1173552da535bfea566f9f5d6f009f14 *data/colonTP.rda
-61d40b7480d7ee79a0c5c7c845dd24c8 *man/TP.idm-package.Rd
-230a1d11527f8c9f9cbc506cc96f2e12 *man/TPidm.Rd
-8d680d14ceb3560a38b7b50dd2dde478 *man/colonTP.Rd
+fd99c1b235af96a89009fe435ad35022 *R/test.nm.R
+020ed21bfd042823c934280564ad5de4 *R/var.AJ.R
+5ea2a433a6ee003343291e6e16dcd747 *R/var.NM.R
+99cc5fd5537a242741fa797ff6db1f40 *build/partial.rdb
+0cbb5cbee64935596b736bbe38ef1c60 *data/colonTP.rda
+3f290b5814bed5cae5713e96bc837f68 *inst/CITATION
+1b691f8285aae4bac744b1c054692f0a *man/TP.idm-package.Rd
+5096eaaaf52efb6e6a9814016bc1e503 *man/TPidm.Rd
+7edfedb9e6b7c8c5ef79fa686043e0ff *man/colonTP.Rd
 6a2b979799934b5e2c4d024e7a7cfd57 *man/plot.TPidm.Rd
 251d82ea5b3153ef6f63649e1e9f580e *man/print.TPidm.Rd
 556686051adf4746e8ee4f12395db551 *man/summary.TPidm.Rd
 dcb0c11d014a8f4f571e7f69d22ceb3e *man/test.nm.Rd
 2221a2d775fba27c2f4e86428e082b2c *src/Cfunctions.c
+9d70ca681fda22d7056bc201031b0107 *src/init.c
diff --git a/NAMESPACE b/NAMESPACE
index 6bdda38..4fc7943 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,4 +1,4 @@
-useDynLib(TP.idm)
+useDynLib(TP.idm, .registration=TRUE)
 importFrom("graphics","layout","layout.show","lines","plot","title","abline","par")
 importFrom("stats","qnorm","var")
 export(plot.TPidm, print.TPidm, summary.TPidm, TPidm, test.nm)
diff --git a/R/TP.idm-internal.R b/R/TP.idm-internal.R
new file mode 100644
index 0000000..113ef12
--- /dev/null
+++ b/R/TP.idm-internal.R
@@ -0,0 +1,126 @@
+.Random.seed <-
+c(403L, 10L, 1317493215L, -1988919883L, -1987163316L, -587781726L, 
+1037527293L, -1379518985L, -1613376738L, 1968459224L, 1626449723L, 
+-1820651303L, -598655896L, -1204420202L, -1887354607L, -1864928573L, 
+-1603040558L, 995097444L, 412053159L, -1297426435L, -2121723404L, 
+1913528026L, 101048357L, 1475996991L, -1133223226L, 2112016208L, 
+-575647885L, 996888977L, 1205371840L, -1508366946L, -1619142807L, 
+-1383967845L, -402269526L, -1655327540L, -954970737L, -881284667L, 
+1731841404L, 137016786L, 742576461L, 79057927L, -719215634L, 
+1846225608L, 1624245707L, -1473601495L, 1212089976L, -1261288538L, 
+933856705L, 1090044499L, 1971030274L, -1442103244L, 29129143L, 
+410054061L, -730749244L, -359019798L, -1520358059L, -40672529L, 
+-382500618L, 62709792L, -891526685L, -1127357823L, -173667088L, 
+443499406L, 2041380089L, -811252085L, 118945466L, 946469308L, 
+-1146868993L, -1460941227L, -205629076L, 1808583106L, -248454627L, 
+1482177175L, -1833167874L, -161386568L, 115819419L, 1702087161L, 
+1112217928L, -1421220746L, 527937969L, 1235902435L, 456326322L, 
+-1856663292L, 1787484743L, -1141826275L, 359230356L, -1923211462L, 
+-2018212987L, 1308112991L, -383808794L, -1229344400L, -1177232749L, 
+90440497L, 1237639648L, 950294142L, -20423415L, 1605021499L, 
+-1876564214L, -68512020L, 1463092079L, -1870187355L, -683776420L, 
+657354290L, 1596493101L, 746470375L, -1947944178L, 666991848L, 
+-1562077461L, -624789623L, 1634681624L, 1841649222L, -2106571103L, 
+-1769993933L, -1562548062L, -239392876L, -334690793L, -999050739L, 
+935020324L, -2010119414L, -1532670347L, -1694388785L, -1077651754L, 
+1730632320L, -882408253L, -85990815L, -1574077744L, -1508156370L, 
+870904281L, 179498475L, 1471998042L, -437442020L, 1228851999L, 
+1226784629L, 300541452L, -1746747422L, 765494461L, 583009079L, 
+-2004112418L, -831887080L, 355507323L, -1343061479L, -1034427736L, 
+-317929386L, 1834337489L, -140722045L, -1869753582L, 736609956L, 
+-1152421273L, -556471363L, -192395980L, -2099528038L, 2015333605L, 
+-1878282625L, 1654175750L, 53640080L, 233375283L, 1545541073L, 
+-1067030656L, 1533049438L, -981655127L, -1649635237L, -51787670L, 
+-885036788L, 1586668367L, 407404549L, 1746067132L, -943032430L, 
+646718605L, 1161200327L, 1765307310L, 1542630920L, -193953269L, 
+-246394135L, 1504468408L, 1407228518L, 1172956929L, 1779683731L, 
+-497194430L, -1413677580L, 326977527L, 335690733L, 2107479172L, 
+1938983594L, 981557269L, 1712687279L, -1214150986L, 978408672L, 
+-1652918237L, -1982863551L, 1338870320L, -954230706L, 645253305L, 
+-1847853877L, 2032189818L, 1227947132L, -1031440193L, 1616498069L, 
+1292668332L, 2140653442L, -2098770595L, 67028055L, -586456258L, 
+659917688L, -617240485L, 392097721L, -174510584L, -1645249354L, 
+2095508465L, -872802781L, -968309134L, 1243421892L, -1950426489L, 
+-267824547L, 148657236L, 2070429306L, -340836923L, -110608609L, 
+-423236954L, -387888080L, -1601245485L, -1825927951L, -354432992L, 
+16232126L, -565392695L, -190693253L, 357432906L, -1551504724L, 
+765648303L, -99373723L, 533674780L, -200934456L, -1036073798L, 
+1748438928L, -1088471572L, -1158553516L, -1852731758L, 83317536L, 
+1685214668L, -80258080L, -697791806L, 754534952L, -29186388L, 
+698886844L, -1318880782L, 362039600L, 748231844L, -1457812808L, 
+-124456102L, 415724944L, 1414480252L, -1932110572L, 476492098L, 
+1739635456L, 1979530172L, 1363927888L, 843634978L, -154616504L, 
+-1893480436L, 2078982236L, -1312256814L, 1704585216L, -715891884L, 
+743664104L, -477667334L, -1266082096L, 503493932L, -228994828L, 
+412850066L, 721435744L, 213429516L, -1766580512L, -2103493278L, 
+-662605144L, -1483929012L, -849860740L, 104725298L, -2058688912L, 
+-110482108L, -1112040488L, 2046574714L, 1026242096L, -1761037636L, 
+344336212L, -758705598L, 1363096480L, 609398972L, 2002692560L, 
+-12667550L, 1905810408L, 1290502092L, -1684941060L, -234758158L, 
+-1945285504L, 1138432244L, 643536008L, -53082822L, 243884880L, 
+538938476L, 1348399252L, 446486866L, 1354024672L, -1134780212L, 
+-107745120L, 68546946L, -1283993368L, -17194772L, -1874608068L, 
+1818239090L, -1995847568L, -1958323676L, -574086600L, -1539452390L, 
+-1595017008L, -1594034628L, 1693701268L, 627614850L, 11988928L, 
+-1982450372L, 1519167120L, 479444002L, -806967288L, 636269708L, 
+1846994972L, 1540539538L, 1420008576L, -1930740716L, -2101267672L, 
+-1445721030L, 1410710736L, 1111198828L, -1220331852L, 1935089682L, 
+1957158880L, 2138149836L, 185890272L, -63759326L, -451799768L, 
+383172876L, 1070603836L, -1043344142L, 757971184L, -740219068L, 
+568885336L, -2143733190L, -295491344L, -1450069316L, 731542036L, 
+-228700990L, 1200424544L, 2108216700L, 221064400L, -2128389086L, 
+645204136L, -1339114100L, -841176644L, 652261426L, -1936221376L, 
+-1363712332L, -370306744L, 596837050L, -1482598640L, -805738260L, 
+-68364204L, 1151758098L, 1311007136L, -627929140L, -639256736L, 
+-1328433982L, 881352232L, -365327956L, -1885015236L, -1904372878L, 
+1461786416L, 1934084132L, 87023160L, -1434111270L, 1266460560L, 
+-664638212L, -906652268L, -1178234046L, -1961874944L, 2013792316L, 
+-1917064624L, 483475618L, 1818923848L, 882328076L, 1934873308L, 
+-1424832430L, 1024232576L, -887640236L, -24420376L, 2021415674L, 
+-445956272L, -1227140052L, -1295199884L, 381572754L, -1684039072L, 
+1743984780L, 2002774368L, -1567110942L, -753466200L, 1075961548L, 
+-1173632516L, 1478046642L, 2030443504L, -1259533884L, -129823400L, 
+1865218426L, -494970448L, 1132839356L, -553324972L, 1106393538L, 
+1196534560L, 588086204L, 1759989840L, -319331358L, -412858520L, 
+-400064692L, -466391300L, -90992910L, -383701760L, 468313204L, 
+764767752L, 415292858L, -1130928688L, 766906348L, 524741524L, 
+610439378L, 956087392L, 669177164L, 49154848L, 603224706L, 133816808L, 
+-1920146452L, -338266564L, -1583178126L, 1868706672L, 80129316L, 
+1835252280L, -1172041318L, -1482455600L, 1972250044L, 1272235668L, 
+-556793214L, 1485444544L, -1979609284L, 594121360L, -1624389854L, 
+2136151560L, 1850819340L, 70523804L, -471931886L, -875252864L, 
+-1789744364L, -288910040L, -2112393158L, -28994608L, -948158762L, 
+932522851L, 26993445L, 841232578L, -1225507408L, 37912025L, 1399814987L, 
+-1171641924L, -100726294L, 939396639L, 2142178473L, 57021502L, 
+1543597356L, 1890277245L, 755832263L, -1930473296L, -750137650L, 
+1145869307L, 730633309L, -581973110L, -144540760L, -1593116399L, 
+1505103875L, -1040411900L, 1336323890L, -1240904313L, -131181615L, 
+1099476310L, 778404660L, 1699577413L, 1834855055L, 860185608L, 
+1326870566L, 1627189459L, -766541259L, -119389806L, -1904155296L, 
+-156283703L, 324825083L, -1334259764L, -1875951014L, 1364148047L, 
+-1156485607L, 369708078L, -1786867140L, -599380051L, -956962729L, 
+566218144L, -1238343938L, -988463669L, 1736813517L, 1473732346L, 
+692344120L, 222215905L, 2067632531L, -707692108L, 1910942402L, 
+1520589527L, 1519061793L, -1615254810L, -1778594844L, 536603413L, 
+-147950785L, -1075439912L, 603214582L, 613084419L, 1065092677L, 
+-1515307422L, 943953360L, -2123963655L, 1752682539L, 377003932L, 
+-1273278518L, -1698764609L, -1528440311L, 1922052894L, -520340788L, 
+-1325223267L, -2015125721L, 1025552720L, -327499794L, 1188503707L, 
+-1828615427L, 832893418L, 1543789192L, -1024447439L, 1685184163L, 
+1387861668L, 1480570706L, -694870361L, -999595279L, -2110313930L, 
+1489643220L, -154940827L, 298731631L, 340863656L, -702220282L, 
+1113568563L, 810613141L, -1905080846L, -1342350784L, -2103037015L, 
+1413459483L, -274651156L, 1411911674L, 1786819631L, -721812423L, 
+-1625622962L, 1345390236L, 1189592077L, 1313505399L, 451381376L, 
+269897438L, 289383851L, -1250789715L, 495536922L, 1487270872L, 
+778612801L, 1856855795L, -816250604L, 1218925730L, 676975031L, 
+1939297665L, 185508742L, -1751715388L, -745080587L, -815294113L, 
+1859359672L, -922831210L, -1482653533L, -1151112603L, 1961820034L, 
+-196832016L, -1872743L, 327943307L, 834473340L, -2138621014L, 
+134478815L, 1335259497L, 1038228862L, -195887764L, 1181132221L, 
+664839815L, -1244142992L, 1656195342L, 2066486843L, 2035664157L, 
+1725485002L, -512501528L, 1916435409L, -590725693L, -1516480828L, 
+-1705910926L, -2019225529L, 435655697L, -1024944746L, 1232234740L, 
+1334076037L, 100223951L, 220299976L, -286751898L, -1345570413L, 
+-1525452043L, 1715432146L, -1324624608L, -1436787447L, -1742640581L, 
+-594690164L, 1492721178L, 484234450L)
diff --git a/R/TPidm.R b/R/TPidm.R
index 8b401a9..c8b60be 100644
--- a/R/TPidm.R
+++ b/R/TPidm.R
@@ -10,7 +10,7 @@ function(data, s, t="last", cov=NULL, CI=TRUE, level=0.95, ci.transformation="li
     } else{
       data<-data[,c("time1","event1","Stime","event",cov)]
       if(length(names(data))!=5) stop("'data' must have 5 variables")
-      if(sum(c("time1","event1","Stime","event",cov)==names(data))!=5) stop("'data' must contain the rigth variables")
+      if(sum(c("time1","event1","Stime","event",cov)==names(data))!=5) stop("'data' must contain the rigth variables")           
       cc<-match(cov,names(data))
       if(!is.factor(data[,cc])) stop("In cov a factor is needed")
     }
@@ -45,7 +45,7 @@ function(data, s, t="last", cov=NULL, CI=TRUE, level=0.95, ci.transformation="li
     # non-absorbing states
     tr.states <- states[!states=="3"]
     
-    if (is.null(cov)){
+      if (is.null(cov)){
       # choose method
       method.type<-c("NM","AJ")
       m<-charmatch(method,method.type,nomatch=0)
@@ -208,7 +208,7 @@ function(data, s, t="last", cov=NULL, CI=TRUE, level=0.95, ci.transformation="li
           if(CI==TRUE){
             times_i<-results.times.NM[[i]]
             probs_i<-results.all.probs.NM[[i]]
-            
+
             
             # var using Titman method:
             variances<- var.NM(data_i, ns, states, tr.states, s, t, probs_i) 
@@ -220,59 +220,59 @@ function(data, s, t="last", cov=NULL, CI=TRUE, level=0.95, ci.transformation="li
             results.ci.NM[[i]]<-ci$CI
           }
         }
-        
+
         
         m.t.nm <-vector('list',n.cat)
         for(i in 1:n.cat){
           m.t.nm[i]<-max(results.times.NM[[i]])
         }
+
+names(results.times.NM) <- paste(rep("t",n.cat), sep="")
+names(results.all.probs.NM) <- names(results.probs.NM) <- names(results.all.est.NM) <- paste(rep("probs",n.cat),sep="")
+names(m.t.nm) <- paste(rep("t",n.cat), sep="")
+
+if(CI==TRUE){
+  names(results.all.ci.NM)<-paste(rep("CI",n.cat), sep="")
+  names(results.ci.NM)<-paste(rep("CI",n.cat), sep="")
+}
         
-        names(results.times.NM) <- paste(rep("t",n.cat), sep="")
-        names(results.all.probs.NM) <- names(results.probs.NM) <- names(results.all.est.NM) <- paste(rep("probs",n.cat),sep="")
-        names(m.t.nm) <- paste(rep("t",n.cat), sep="")
-        
-        if(CI==TRUE){
-          names(results.all.ci.NM)<-paste(rep("CI",n.cat), sep="")
-          names(results.ci.NM)<-paste(rep("CI",n.cat), sep="")
-        }
-        
-        if (s==0){
-          p.trans<-c("1 1", "1 2", "1 3")
-        }else{
-          p.trans<-c("1 1", "1 2", "1 3", "2 2", "2 3")
-        }
-        
-        if(CI==TRUE){
-          # results:
-          res <- list(
-            # states information:
-            cov=cov, names.cov=names.cov,
-            method=method,s=s,
-            t=m.t.nm,
-            states=states, ns=ns, tr.states=tr.states, 
-            ci.transformation=ci.transformation,
-            # event times:
-            times=results.times.NM,
-            # confidence intervals:
-            probs=results.ci.NM, all.probs=results.all.ci.NM,
-            # posible transitions:
-            p.trans=p.trans,CI=CI)
-        }else{
-          # results:
-          res <- list(
-            # states information:
-            cov=cov, names.cov=names.cov,
-            method=method,s=s,
-            t=m.t.nm,
-            states=states, ns=ns, tr.states=tr.states,
-            ci.transformation=ci.transformation,
-            # event times:
-            times=results.times.NM,
-            # occupation or transition probabilities:
-            probs=results.probs.NM, all.probs=results.all.est.NM,
-            # posible transitions:
-            p.trans=p.trans,CI=CI)
-        }
+if (s==0){
+  p.trans<-c("1 1", "1 2", "1 3")
+}else{
+  p.trans<-c("1 1", "1 2", "1 3", "2 2", "2 3")
+}
+
+if(CI==TRUE){
+  # results:
+  res <- list(
+    # states information:
+    cov=cov, names.cov=names.cov,
+    method=method,s=s,
+    t=m.t.nm,
+    states=states, ns=ns, tr.states=tr.states, 
+    ci.transformation=ci.transformation,
+    # event times:
+    times=results.times.NM,
+    # confidence intervals:
+    probs=results.ci.NM, all.probs=results.all.ci.NM,
+    # posible transitions:
+    p.trans=p.trans,CI=CI)
+}else{
+  # results:
+  res <- list(
+    # states information:
+    cov=cov, names.cov=names.cov,
+    method=method,s=s,
+    t=m.t.nm,
+    states=states, ns=ns, tr.states=tr.states,
+    ci.transformation=ci.transformation,
+    # event times:
+    times=results.times.NM,
+    # occupation or transition probabilities:
+    probs=results.probs.NM, all.probs=results.all.est.NM,
+    # posible transitions:
+    p.trans=p.trans,CI=CI)
+}
       }
       if (m==2){
         
@@ -344,62 +344,62 @@ function(data, s, t="last", cov=NULL, CI=TRUE, level=0.95, ci.transformation="li
             results.ci.AJ[[ii]]<-ci$CI
           }
         }
-        
+              
         m.t.aj <-vector('list',n.cat)
         for(i in 1:n.cat){
           m.t.aj[i]<-max(results.times.AJ[[i]])
         }
-        
-        names(results.times.AJ) <- paste(rep("t",n.cat), sep="")
-        names(results.probs.AJ) <- names(results.all.est.AJ) <- paste(rep("probs",n.cat), sep="")
-        names(m.t.aj) <- paste(rep("t",n.cat), sep="")
-        
-        if(CI==TRUE){
-          names(results.all.ci.AJ)<-paste(rep("CI",n.cat), sep="")
-          names(results.ci.AJ)<-paste(rep("CI",n.cat), sep="")
-        }
-        
-        if (s==0){
-          p.trans<-c("1 1", "1 2", "1 3")
-        }else{
-          p.trans<-c("1 1", "1 2", "1 3", "2 2", "2 3")
-        } 
-        
-        
-        if(CI==TRUE){
-          # results:
-          res <- list(
-            # states information:
-            cov=cov, names.cov=names.cov,
-            method=method,s=s,
-            t=m.t.aj,
-            states=states, ns=ns, tr.states=tr.states, 
-            ci.transformation=ci.transformation,
-            # event times:
-            times=results.times.AJ,
-            # confidence intervals:
-            probs=results.ci.AJ, all.probs=results.all.ci.AJ,
-            # posible transitions:
-            p.trans=p.trans,CI=CI)
-        }else{
-          # results:
-          res <- list(
-            # states information:
-            cov=cov, names.cov=names.cov,
-            method=method,s=s,
-            t=m.t.aj,
-            states=states, ns=ns, tr.states=tr.states,
-            ci.transformation=ci.transformation,
-            # event times:
-            times=results.times.AJ,
-            # occupation or transition probabilities:
-            probs=results.probs.AJ, all.probs=results.all.est.AJ,
-            # posible transitions:
-            p.trans=p.trans,CI=CI)
-        }      
+
+names(results.times.AJ) <- paste(rep("t",n.cat), sep="")
+names(results.probs.AJ) <- names(results.all.est.AJ) <- paste(rep("probs",n.cat), sep="")
+names(m.t.aj) <- paste(rep("t",n.cat), sep="")
+
+if(CI==TRUE){
+  names(results.all.ci.AJ)<-paste(rep("CI",n.cat), sep="")
+  names(results.ci.AJ)<-paste(rep("CI",n.cat), sep="")
+}
+
+if (s==0){
+  p.trans<-c("1 1", "1 2", "1 3")
+}else{
+  p.trans<-c("1 1", "1 2", "1 3", "2 2", "2 3")
+} 
+
+
+if(CI==TRUE){
+  # results:
+  res <- list(
+    # states information:
+    cov=cov, names.cov=names.cov,
+    method=method,s=s,
+    t=m.t.aj,
+    states=states, ns=ns, tr.states=tr.states, 
+    ci.transformation=ci.transformation,
+    # event times:
+    times=results.times.AJ,
+    # confidence intervals:
+    probs=results.ci.AJ, all.probs=results.all.ci.AJ,
+    # posible transitions:
+    p.trans=p.trans,CI=CI)
+}else{
+  # results:
+  res <- list(
+    # states information:
+    cov=cov, names.cov=names.cov,
+    method=method,s=s,
+    t=m.t.aj,
+    states=states, ns=ns, tr.states=tr.states,
+    ci.transformation=ci.transformation,
+    # event times:
+    times=results.times.AJ,
+    # occupation or transition probabilities:
+    probs=results.probs.AJ, all.probs=results.all.est.AJ,
+    # posible transitions:
+    p.trans=p.trans,CI=CI)
+}      
       }
-    }
-    
+      }
+
     
     res$call<-match.call()
     class(res) = "TPidm"
diff --git a/R/ci.AJ.R b/R/ci.AJ.R
index b95c20c..049c3d5 100644
--- a/R/ci.AJ.R
+++ b/R/ci.AJ.R
@@ -1,106 +1,106 @@
 ci.AJ <-
 function(s,t,level,ci.transformation="linear",dNs.id_tr,TP.AJs,cov.AJs,e.times.id_tr){
+  
+  if(s==0){
+    p.transitions <- c("1 1", "1 2", "1 3")
     
-    if(s==0){
-      p.transitions <- c("1 1", "1 2", "1 3")
-      
-      alpha<-qnorm(level+(1-level)/2)
-      
-      CI.tp<-array(0,dim=c(nrow(dNs.id_tr),4,length(p.transitions)),
-                   dimnames=list(rows=e.times.id_tr,cols=c("probs","lower","upper","variance"),trans=p.transitions))
-      
-      # different transformations to built CI
-      ci.transformation <- match.arg(ci.transformation, c("linear", "log", "cloglog", "log-log"))
-      
-      for (j in 1:length(p.transitions)) { ## loop through possible transitions
-        
-        idx <- unlist(strsplit(p.transitions[j], " "))
-        CI.tp[ , 1, j] <- P <- TP.AJs[idx[1], idx[2] , ]
-        CI.tp[ , 4, j] <- var <- cov.AJs[p.transitions[j], p.transitions[j], ]
-        
-        
-        switch(ci.transformation[1],
-               "linear" = {
-                 CI.tp[ , 2, j] <- P - alpha * sqrt(var)
-                 CI.tp[ , 3, j] <- P + alpha * sqrt(var)},
-               "log" = {
-                 CI.tp[ , 2, j] <- exp(log(P) - alpha * sqrt(var) / P)
-                 CI.tp[ , 3, j] <- exp(log(P) + alpha * sqrt(var) / P)},
-               "cloglog" = {
-                 CI.tp[ , 2, j] <- 1 - (1 - P)^(exp(alpha * (sqrt(var) /
-                                                               ((1 - P) * log(1 - P)))))
-                 CI.tp[ , 3, j] <- 1 - (1 - P)^(exp(-alpha * (sqrt(var) /
-                                                                ((1 - P) * log(1 - P)))))},
-               "log-log" = {
-                 CI.tp[ , 2, j] <- P^(exp(-alpha * (sqrt(var) / (P * log(P)))))
-                 CI.tp[ , 3, j] <- P^(exp(alpha * (sqrt(var) / (P * log(P)))))})
-        
-        CI.tp[ , 2, j] <- pmax(CI.tp[ , 2, j], 0)
-        CI.tp[ , 3, j] <- pmin(CI.tp[ , 3, j], 1)
-        
-      } ## end j loop
+    alpha<-qnorm(level+(1-level)/2)
+    
+    CI.tp<-array(0,dim=c(nrow(dNs.id_tr),4,length(p.transitions)),
+                 dimnames=list(rows=e.times.id_tr,cols=c("probs","lower","upper","variance"),trans=p.transitions))
+    
+    # different transformations to built CI
+    ci.transformation <- match.arg(ci.transformation, c("linear", "log", "cloglog", "log-log"))
+    
+    for (j in 1:length(p.transitions)) { ## loop through possible transitions
       
-      CI.t<- matrix(0, nrow=length(p.transitions), ncol=4)
-      colnames(CI.t) <- c("probs","lower","upper","variance")
-      rownames(CI.t) <- p.transitions
+      idx <- unlist(strsplit(p.transitions[j], " "))
+      CI.tp[ , 1, j] <- P <- TP.AJs[idx[1], idx[2] , ]
+      CI.tp[ , 4, j] <- var <- cov.AJs[p.transitions[j], p.transitions[j], ]
       
-      for(j in 1:length(p.transitions)){
-        CI.t[j,]<-CI.tp[nrow(CI.tp[,,j]),,j] 
-      }
       
-      #CI.tp <- round(CI.tp,7)
-      #CI.t <- round(CI.t,7)
-    }else{
-      p.transitions <- c("1 1", "1 2", "1 3", "2 2", "2 3")
+      switch(ci.transformation[1],
+             "linear" = {
+               CI.tp[ , 2, j] <- P - alpha * sqrt(var)
+               CI.tp[ , 3, j] <- P + alpha * sqrt(var)},
+             "log" = {
+               CI.tp[ , 2, j] <- exp(log(P) - alpha * sqrt(var) / P)
+               CI.tp[ , 3, j] <- exp(log(P) + alpha * sqrt(var) / P)},
+             "cloglog" = {
+               CI.tp[ , 2, j] <- 1 - (1 - P)^(exp(alpha * (sqrt(var) /
+                                                             ((1 - P) * log(1 - P)))))
+               CI.tp[ , 3, j] <- 1 - (1 - P)^(exp(-alpha * (sqrt(var) /
+                                                              ((1 - P) * log(1 - P)))))},
+             "log-log" = {
+               CI.tp[ , 2, j] <- P^(exp(-alpha * (sqrt(var) / (P * log(P)))))
+               CI.tp[ , 3, j] <- P^(exp(alpha * (sqrt(var) / (P * log(P)))))})
       
-      alpha<-qnorm(level+(1-level)/2)
+      CI.tp[ , 2, j] <- pmax(CI.tp[ , 2, j], 0)
+      CI.tp[ , 3, j] <- pmin(CI.tp[ , 3, j], 1)
       
-      CI.tp<-array(0,dim=c(nrow(dNs.id_tr),4,length(p.transitions)),
-                   dimnames=list(rows=e.times.id_tr,cols=c("probs","lower","upper","variance"),trans=p.transitions))
+    } ## end j loop
+    
+    CI.t<- matrix(0, nrow=length(p.transitions), ncol=4)
+    colnames(CI.t) <- c("probs","lower","upper","variance")
+    rownames(CI.t) <- p.transitions
+    
+    for(j in 1:length(p.transitions)){
+      CI.t[j,]<-CI.tp[nrow(CI.tp[,,j]),,j] 
+    }
+    
+    #CI.tp <- round(CI.tp,7)
+    #CI.t <- round(CI.t,7)
+  }else{
+    p.transitions <- c("1 1", "1 2", "1 3", "2 2", "2 3")
+    
+    alpha<-qnorm(level+(1-level)/2)
+    
+    CI.tp<-array(0,dim=c(nrow(dNs.id_tr),4,length(p.transitions)),
+                 dimnames=list(rows=e.times.id_tr,cols=c("probs","lower","upper","variance"),trans=p.transitions))
+    
+    # different transformations to built CI
+    ci.transformation <- match.arg(ci.transformation, c("linear", "log", "cloglog", "log-log"))
+    
+    for (j in 1:length(p.transitions)) { ## loop through possible transitions
       
-      # different transformations to built CI
-      ci.transformation <- match.arg(ci.transformation, c("linear", "log", "cloglog", "log-log"))
+      idx <- unlist(strsplit(p.transitions[j], " "))
+      CI.tp[ , 1, j] <- P <- TP.AJs[idx[1], idx[2] , ]
+      CI.tp[ , 4, j] <- var <- cov.AJs[p.transitions[j], p.transitions[j], ]
       
-      for (j in 1:length(p.transitions)) { ## loop through possible transitions
-        
-        idx <- unlist(strsplit(p.transitions[j], " "))
-        CI.tp[ , 1, j] <- P <- TP.AJs[idx[1], idx[2] , ]
-        CI.tp[ , 4, j] <- var <- cov.AJs[p.transitions[j], p.transitions[j], ]
-        
-        
-        switch(ci.transformation[1],
-               "linear" = {
-                 CI.tp[ , 2, j] <- P - alpha * sqrt(var)
-                 CI.tp[ , 3, j] <- P + alpha * sqrt(var)},
-               "log" = {
-                 CI.tp[ , 2, j] <- exp(log(P) - alpha * sqrt(var) / P)
-                 CI.tp[ , 3, j] <- exp(log(P) + alpha * sqrt(var) / P)},
-               "cloglog" = {
-                 CI.tp[ , 2, j] <- 1 - (1 - P)^(exp(alpha * (sqrt(var) /
-                                                               ((1 - P) * log(1 - P)))))
-                 CI.tp[ , 3, j] <- 1 - (1 - P)^(exp(-alpha * (sqrt(var) /
-                                                                ((1 - P) * log(1 - P)))))},
-               "log-log" = {
-                 CI.tp[ , 2, j] <- P^(exp(-alpha * (sqrt(var) / (P * log(P)))))
-                 CI.tp[ , 3, j] <- P^(exp(alpha * (sqrt(var) / (P * log(P)))))})
-        
-        CI.tp[ , 2, j] <- pmax(CI.tp[ , 2, j], 0)
-        CI.tp[ , 3, j] <- pmin(CI.tp[ , 3, j], 1)
-        
-      } ## end j loop
       
-      CI.t<- matrix(0, nrow=length(p.transitions), ncol=4)
-      colnames(CI.t) <- c("probs","lower","upper","variance")
-      rownames(CI.t) <- p.transitions
+      switch(ci.transformation[1],
+             "linear" = {
+               CI.tp[ , 2, j] <- P - alpha * sqrt(var)
+               CI.tp[ , 3, j] <- P + alpha * sqrt(var)},
+             "log" = {
+               CI.tp[ , 2, j] <- exp(log(P) - alpha * sqrt(var) / P)
+               CI.tp[ , 3, j] <- exp(log(P) + alpha * sqrt(var) / P)},
+             "cloglog" = {
+               CI.tp[ , 2, j] <- 1 - (1 - P)^(exp(alpha * (sqrt(var) /
+                                                             ((1 - P) * log(1 - P)))))
+               CI.tp[ , 3, j] <- 1 - (1 - P)^(exp(-alpha * (sqrt(var) /
+                                                              ((1 - P) * log(1 - P)))))},
+             "log-log" = {
+               CI.tp[ , 2, j] <- P^(exp(-alpha * (sqrt(var) / (P * log(P)))))
+               CI.tp[ , 3, j] <- P^(exp(alpha * (sqrt(var) / (P * log(P)))))})
       
-      for(j in 1:length(p.transitions)){
-        CI.t[j,]<-CI.tp[nrow(CI.tp[,,j]),,j] 
-      }
+      CI.tp[ , 2, j] <- pmax(CI.tp[ , 2, j], 0)
+      CI.tp[ , 3, j] <- pmin(CI.tp[ , 3, j], 1)
       
-      #CI.tp <- round(CI.tp,7)
-      #CI.t <-round(CI.t,7)
-    }
+    } ## end j loop
     
-    return(list(CI=CI.t,all.CI=CI.tp))
+    CI.t<- matrix(0, nrow=length(p.transitions), ncol=4)
+    colnames(CI.t) <- c("probs","lower","upper","variance")
+    rownames(CI.t) <- p.transitions
+    
+    for(j in 1:length(p.transitions)){
+      CI.t[j,]<-CI.tp[nrow(CI.tp[,,j]),,j] 
+    }
     
+    #CI.tp <- round(CI.tp,7)
+    #CI.t <-round(CI.t,7)
   }
+  
+  return(list(CI=CI.t,all.CI=CI.tp))
+  
+}
diff --git a/R/fun.AJ.R b/R/fun.AJ.R
index 1ec1f0f..579e28a 100644
--- a/R/fun.AJ.R
+++ b/R/fun.AJ.R
@@ -1,99 +1,99 @@
 fun.AJ <-
 function(ns,states,dNs,Ys,sum_dNs,s,t, event.times,initial.probs){
+  
+  if (t=="last") t <- event.times[length(event.times)]
+  
+  id_tr <- which(s<event.times& event.times<=t) ## location of those (s, t]
+  
+  l.id_tr <- length(id_tr)
+  
+  # to use in variance function:
+  e.times.id_tr<-event.times[id_tr]
+  
+  ## Occupation Probability matrix
+  OP <- matrix(NA, nrow=l.id_tr, ncol=length(states))
+  rownames(OP) <- rownames(dNs)[id_tr]; colnames(OP) <- states
+  all.dA <- all.I.dA <- array(dim=c(length(states), length(states), l.id_tr),
+                              dimnames=list(rows=states,cols=states, time=rownames(dNs)[id_tr]))
+  
+  ## Transition Probability matrix
+  cum.prod <- diag(ns)
+  rownames(cum.prod) <- states
+  
+  TP.AJs <- array(dim=c(ns, ns, l.id_tr),dimnames=list(rows=states, cols=states,
+                                                       dim=rownames(dNs)[id_tr]))
+  
+  dNs.id_tr<-dNs[id_tr,]
+  Ys.id_tr<-Ys[id_tr,]
+  sum_dNs.id_tr<-sum_dNs[id_tr,]
+  
+  
+  for(i in 1:l.id_tr){
     
-    if (t=="last") t <- event.times[length(event.times)]
-    
-    id_tr <- which(s<event.times& event.times<=t) ## location of those (s, t]
-    
-    l.id_tr <- length(id_tr)
-    
-    # to use in variance function:
-    e.times.id_tr<-event.times[id_tr]
-    
-    ## Occupation Probability matrix
-    OP <- matrix(NA, nrow=l.id_tr, ncol=length(states))
-    rownames(OP) <- rownames(dNs)[id_tr]; colnames(OP) <- states
-    all.dA <- all.I.dA <- array(dim=c(length(states), length(states), l.id_tr),
-                                dimnames=list(rows=states,cols=states, time=rownames(dNs)[id_tr]))
-    
-    ## Transition Probability matrix
-    cum.prod <- diag(ns)
-    rownames(cum.prod) <- states
-    
-    TP.AJs <- array(dim=c(ns, ns, l.id_tr),dimnames=list(rows=states, cols=states,
-                                                         dim=rownames(dNs)[id_tr]))
-    
-    dNs.id_tr<-dNs[id_tr,]
-    Ys.id_tr<-Ys[id_tr,]
-    sum_dNs.id_tr<-sum_dNs[id_tr,]
-    
-    
-    for(i in 1:l.id_tr){
-      
-      I.dA <- diag(ns) ## creates trans matrix for current time
-      
-      dA <- matrix(0, nrow=ns, ncol=ns)
-      
-      colnames(I.dA) <- rownames(I.dA) <- colnames(dA) <- rownames(dA) <- states
-      
-      i_tr <- which(dNs.id_tr[i, , drop=FALSE]>0)  ## indicator for kind of transition at time i
-      dNs.event.names <- colnames(dNs.id_tr)[i_tr] ## gets names of transitions (ie:  dN##)
-      split_dNs.event <- strsplit(dNs.event.names, " ")    ## splits title of dN##
-      st.start <- sapply(split_dNs.event, function(x) x[2])
-      st.end <- sapply(split_dNs.event, function(x) x[3])  ## start & stop states as character strings
-      i_tr.s <- matrix(as.character(c(st.start, st.end)), ncol=2)
-      i_tr.s2 <- matrix(as.character(c(st.start, st.start)), ncol=2)
-      
-      dA[i_tr.s] <- dNs.id_tr[i, i_tr]/Ys.id_tr[i, paste("Y", st.start)]
-      if (length(i_tr)==1) {
-        dA[st.start, st.start] <- -dNs.id_tr[i, i_tr]/Ys.id_tr[i, paste("Y", st.start)]
-      } else {
-        dA[i_tr.s2] <- -rowSums(dA[st.start, ])
-      }
-      
-      I.dA <- I.dA + dA ## I+dA (transition) matrix
-      
-      all.dA[, , i] <- dA     ## stores all dA matrices
-      all.I.dA[, , i] <- I.dA ## array for storing all tran matrices
-      
-      cum.prod <- cum.prod %*% I.dA
-      TP.AJs[,,i] <- cum.prod
-    } ## end loop
-    
-    if (s==0){    
-      OP[i, ] <- initial.probs%*%TP.AJs[, , i] ## state occupation probabilities 
-      op<-OP[i,]   
-      
-      p.transitions<-c("1 1", "1 2", "1 3")
-      
-      all.est<-array(0,dim=c(l.id_tr,1,length(p.transitions)),
-                     dimnames=list(rows=rownames(dNs)[id_tr],cols="probs",trans=p.transitions))
-      
-      for(j in 1:length(p.transitions)){
-        idx<-unlist(strsplit(p.transitions[j]," "))
-        all.est[,1,j]<-TP.AJs[idx[1],idx[2],]
-      }
-      
-      
-      res <- list(probs=op,all.est=all.est,TP.AJs=TP.AJs,all.I.dA=all.I.dA, dNs.id_tr=dNs.id_tr, 
-                  Ys.id_tr=Ys.id_tr, sum_dNs.id_tr=sum_dNs.id_tr, e.times.id_tr=e.times.id_tr,p.trans=p.transitions,t=t)
-      return(res)
-    } else{
-      p.transitions<-c("1 1", "1 2", "1 3", "2 2", "2 3") 
-      
-      all.est<-array(0,dim=c(l.id_tr,1,length(p.transitions)),
-                     dimnames=list(rows=rownames(dNs)[id_tr],cols="probs",trans=p.transitions))
-      
-      for(j in 1:length(p.transitions)){
-        idx<-unlist(strsplit(p.transitions[j]," "))
-        all.est[,1,j]<-TP.AJs[idx[1],idx[2],]
-      }
-      
-      
-      res <- list(probs=cum.prod,all.est=all.est,TP.AJs=TP.AJs,all.I.dA=all.I.dA, dNs.id_tr=dNs.id_tr, 
-                  Ys.id_tr=Ys.id_tr, sum_dNs.id_tr=sum_dNs.id_tr, e.times.id_tr=e.times.id_tr,p.trans=p.transitions,t=t)
-      return(res)
+    I.dA <- diag(ns) ## creates trans matrix for current time
+    
+    dA <- matrix(0, nrow=ns, ncol=ns)
+    
+    colnames(I.dA) <- rownames(I.dA) <- colnames(dA) <- rownames(dA) <- states
+    
+    i_tr <- which(dNs.id_tr[i, , drop=FALSE]>0)  ## indicator for kind of transition at time i
+    dNs.event.names <- colnames(dNs.id_tr)[i_tr] ## gets names of transitions (ie:  dN##)
+    split_dNs.event <- strsplit(dNs.event.names, " ")    ## splits title of dN##
+    st.start <- sapply(split_dNs.event, function(x) x[2])
+    st.end <- sapply(split_dNs.event, function(x) x[3])  ## start & stop states as character strings
+    i_tr.s <- matrix(as.character(c(st.start, st.end)), ncol=2)
+    i_tr.s2 <- matrix(as.character(c(st.start, st.start)), ncol=2)
+    
+    dA[i_tr.s] <- dNs.id_tr[i, i_tr]/Ys.id_tr[i, paste("Y", st.start)]
+    if (length(i_tr)==1) {
+      dA[st.start, st.start] <- -dNs.id_tr[i, i_tr]/Ys.id_tr[i, paste("Y", st.start)]
+    } else {
+      dA[i_tr.s2] <- -rowSums(dA[st.start, ])
+    }
+    
+    I.dA <- I.dA + dA ## I+dA (transition) matrix
+    
+    all.dA[, , i] <- dA     ## stores all dA matrices
+    all.I.dA[, , i] <- I.dA ## array for storing all tran matrices
+    
+    cum.prod <- cum.prod %*% I.dA
+    TP.AJs[,,i] <- cum.prod
+  } ## end loop
+  
+  if (s==0){    
+    OP[i, ] <- initial.probs%*%TP.AJs[, , i] ## state occupation probabilities 
+    op<-OP[i,]   
+    
+    p.transitions<-c("1 1", "1 2", "1 3")
+    
+    all.est<-array(0,dim=c(l.id_tr,1,length(p.transitions)),
+                   dimnames=list(rows=rownames(dNs)[id_tr],cols="probs",trans=p.transitions))
+    
+    for(j in 1:length(p.transitions)){
+      idx<-unlist(strsplit(p.transitions[j]," "))
+      all.est[,1,j]<-TP.AJs[idx[1],idx[2],]
+    }
+    
+    
+    res <- list(probs=op,all.est=all.est,TP.AJs=TP.AJs,all.I.dA=all.I.dA, dNs.id_tr=dNs.id_tr, 
+                Ys.id_tr=Ys.id_tr, sum_dNs.id_tr=sum_dNs.id_tr, e.times.id_tr=e.times.id_tr,p.trans=p.transitions,t=t)
+    return(res)
+  } else{
+    p.transitions<-c("1 1", "1 2", "1 3", "2 2", "2 3") 
+    
+    all.est<-array(0,dim=c(l.id_tr,1,length(p.transitions)),
+                   dimnames=list(rows=rownames(dNs)[id_tr],cols="probs",trans=p.transitions))
+    
+    for(j in 1:length(p.transitions)){
+      idx<-unlist(strsplit(p.transitions[j]," "))
+      all.est[,1,j]<-TP.AJs[idx[1],idx[2],]
     }
     
     
+    res <- list(probs=cum.prod,all.est=all.est,TP.AJs=TP.AJs,all.I.dA=all.I.dA, dNs.id_tr=dNs.id_tr, 
+                Ys.id_tr=Ys.id_tr, sum_dNs.id_tr=sum_dNs.id_tr, e.times.id_tr=e.times.id_tr,p.trans=p.transitions,t=t)
+    return(res)
   }
+  
+  
+}
diff --git a/R/plot.TPidm.R b/R/plot.TPidm.R
index dba3fbc..2b4ee27 100644
--- a/R/plot.TPidm.R
+++ b/R/plot.TPidm.R
@@ -11,111 +11,111 @@ function(x,chosen.tr="ALL",col="black", ...){
     
     
     if(is.null(x$cov)){
-      
+    
       if(is.null(col)) col<-c("black","red")
       if(length(col)<2) col<-rep(col,length.out=2)
       
-      # CI==TRUE
-      if(CI==TRUE){
+    # CI==TRUE
+    if(CI==TRUE){
+
+      all.probs<-x$all.probs[,,chosen.tr]
+      
+      # option 1: individual curves with CI
+      if(length(chosen.tr)==1){
+        tit<-paste("p",chosen.tr)
+        timemax<-max(times)*1.04
+        plot(times,all.probs[,"probs"],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),col=col[1],xlim=c(0,timemax),cex.axis=0.7,cex.lab=0.7,main=tit, ...)
+        lines(times,all.probs[,"lower"],type="s",lty=3,col=col[2])
+        lines(times,all.probs[,"upper"],type="s",lty=3,col=col[2])
         
-        all.probs<-x$all.probs[,,chosen.tr]
+      }else if(length(chosen.tr)>1 & length(chosen.tr)<=3){
+        N<-seq(1,length(chosen.tr),1)
+        tpg<-layout(matrix(N,1,length(chosen.tr),byrow=TRUE))
+        layout.show(tpg)
         
-        # option 1: individual curves with CI
-        if(length(chosen.tr)==1){
-          tit<-paste("p",chosen.tr)
+        for(i in 1:length(chosen.tr)){
+          tit<-paste("p",chosen.tr[[i]])
           timemax<-max(times)*1.04
-          plot(times,all.probs[,"probs"],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),col=col[1],xlim=c(0,timemax),cex.axis=0.7,cex.lab=0.7,main=tit, ...)
-          lines(times,all.probs[,"lower"],type="s",lty=3,col=col[2])
-          lines(times,all.probs[,"upper"],type="s",lty=3,col=col[2])
-          
-        }else if(length(chosen.tr)>1 & length(chosen.tr)<=3){
-          N<-seq(1,length(chosen.tr),1)
-          tpg<-layout(matrix(N,1,length(chosen.tr),byrow=TRUE))
-          layout.show(tpg)
-          
-          for(i in 1:length(chosen.tr)){
-            tit<-paste("p",chosen.tr[[i]])
-            timemax<-max(times)*1.04
-            plot(times,all.probs[,"probs",chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
-            lines(times,all.probs[,"lower",chosen.tr[[i]]],type="s",lty=3,col=col[2])
-            lines(times,all.probs[,"upper",chosen.tr[[i]]],type="s",lty=3,col=col[2]) 
-          }
-        } else {
-          
-          n.cols<-3
-          n.rows<-ceiling(length(chosen.tr)/n.cols)
-          N<-seq(1,length(chosen.tr),1)
-          if(length(N)<n.cols*n.rows){
-            m<-n.cols*n.rows-length(N)
-            N<-c(N,rep(0,m))
-          }
-          tpg<-layout(matrix(N,n.rows,n.cols,byrow=TRUE))
-          layout.show(tpg)
-          
-          for(i in 1:length(chosen.tr)){
-            tit<-paste("p",chosen.tr[[i]])
-            timemax<-max(times)*1.04
-            plot(times,all.probs[,"probs",chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",col=col[1],ylim=c(0,1),xlim=c(0,timemax),cex.axis=0.7,cex.lab=0.7,main=tit, ...)
-            lines(times,all.probs[,"lower",chosen.tr[[i]]],type="s",lty=3,col=col[2])
-            lines(times,all.probs[,"upper",chosen.tr[[i]]],type="s",lty=3,col=col[2]) 
-          }
+          plot(times,all.probs[,"probs",chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
+          lines(times,all.probs[,"lower",chosen.tr[[i]]],type="s",lty=3,col=col[2])
+          lines(times,all.probs[,"upper",chosen.tr[[i]]],type="s",lty=3,col=col[2]) 
         }
+      } else {
         
-        if(s==0){
-          title("State Occupation Probabilities", outer=TRUE,  line=-1)
-        } else{
-          title("Transition Probabilities", outer=TRUE,  line=-1)
+        n.cols<-3
+        n.rows<-ceiling(length(chosen.tr)/n.cols)
+        N<-seq(1,length(chosen.tr),1)
+        if(length(N)<n.cols*n.rows){
+          m<-n.cols*n.rows-length(N)
+          N<-c(N,rep(0,m))
         }
+        tpg<-layout(matrix(N,n.rows,n.cols,byrow=TRUE))
+        layout.show(tpg)
         
-      }else{
-        
-        
-        all.probs<-x$all.probs[,,chosen.tr]
+        for(i in 1:length(chosen.tr)){
+          tit<-paste("p",chosen.tr[[i]])
+          timemax<-max(times)*1.04
+          plot(times,all.probs[,"probs",chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",col=col[1],ylim=c(0,1),xlim=c(0,timemax),cex.axis=0.7,cex.lab=0.7,main=tit, ...)
+          lines(times,all.probs[,"lower",chosen.tr[[i]]],type="s",lty=3,col=col[2])
+          lines(times,all.probs[,"upper",chosen.tr[[i]]],type="s",lty=3,col=col[2]) 
+        }
+      }
+      
+      if(s==0){
+        title("State Occupation Probabilities", outer=TRUE,  line=-1)
+      } else{
+        title("Transition Probabilities", outer=TRUE,  line=-1)
+      }
+      
+    }else{
+      
+      
+      all.probs<-x$all.probs[,,chosen.tr]
+      
+      # option 2: individual curves without CI
+      
+      if(length(chosen.tr)==1){
+        tit<-paste("p",chosen.tr)
+        timemax<-max(times)*1.04
+        plot(times,all.probs,type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
         
-        # option 2: individual curves without CI
+      }else if(length(chosen.tr)>1 & length(chosen.tr)<=3){
+        N<-seq(1,length(chosen.tr),1)
+        tpg<-layout(matrix(N,1,length(chosen.tr),byrow=TRUE))
+        layout.show(tpg)
         
-        if(length(chosen.tr)==1){
-          tit<-paste("p",chosen.tr)
+        for(i in 1:length(chosen.tr)){
+          tit<-paste("p",chosen.tr[[i]])
           timemax<-max(times)*1.04
-          plot(times,all.probs,type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
-          
-        }else if(length(chosen.tr)>1 & length(chosen.tr)<=3){
-          N<-seq(1,length(chosen.tr),1)
-          tpg<-layout(matrix(N,1,length(chosen.tr),byrow=TRUE))
-          layout.show(tpg)
+          plot(times,all.probs[,chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
           
-          for(i in 1:length(chosen.tr)){
-            tit<-paste("p",chosen.tr[[i]])
-            timemax<-max(times)*1.04
-            plot(times,all.probs[,chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
-            
-          }
-        } else {
-          
-          n.cols<-3
-          n.rows<-ceiling(length(chosen.tr)/n.cols)
-          N<-seq(1,length(chosen.tr),1)
-          if(length(N)<n.cols*n.rows){
-            m<-n.cols*n.rows-length(N)
-            N<-c(N,rep(0,m))
-          }
-          tpg<-layout(matrix(N,n.rows,n.cols,byrow=TRUE))
-          layout.show(tpg)
-          
-          for(i in 1:length(chosen.tr)){
-            tit<-paste("p",chosen.tr[[i]])
-            timemax<-max(times)*1.04
-            plot(times,all.probs[,chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
-            
-          }
         }
-        if(s==0){
-          title("State Occupation Probabilities", outer=TRUE,  line=-1)
-        } else{
-          title("Transition Probabilities", outer=TRUE,  line=-1)
+      } else {
+        
+        n.cols<-3
+        n.rows<-ceiling(length(chosen.tr)/n.cols)
+        N<-seq(1,length(chosen.tr),1)
+        if(length(N)<n.cols*n.rows){
+          m<-n.cols*n.rows-length(N)
+          N<-c(N,rep(0,m))
         }
+        tpg<-layout(matrix(N,n.rows,n.cols,byrow=TRUE))
+        layout.show(tpg)
         
-      } # if cov=="null"
+        for(i in 1:length(chosen.tr)){
+          tit<-paste("p",chosen.tr[[i]])
+          timemax<-max(times)*1.04
+          plot(times,all.probs[,chosen.tr[[i]]],type="s",xlab="Time",ylab="Probability",ylim=c(0,1),xlim=c(0,timemax),col=col[1],cex.axis=0.7,cex.lab=0.7,main=tit, ...)
+          
+        }
+      }
+      if(s==0){
+        title("State Occupation Probabilities", outer=TRUE,  line=-1)
+      } else{
+        title("Transition Probabilities", outer=TRUE,  line=-1)
+      }
+      
+    } # if cov=="null"
     } else {
       
       maxims<-c(NA)
@@ -144,7 +144,7 @@ function(x,chosen.tr="ALL",col="black", ...){
             lines(x$times[j]$t,x$all.probs[j]$CI[,"lower",chosen.tr[[1]]],type="s",lty=3,col=col[j])
             lines(x$times[j]$t,x$all.probs[j]$CI[,"upper",chosen.tr[[1]]],type="s",lty=3,col=col[j])
           }
-          
+
           
         }else if(length(chosen.tr)>1 & length(chosen.tr)<=3){
           N<-seq(1,length(chosen.tr),1)
diff --git a/R/prep.data.AJ.R b/R/prep.data.AJ.R
index 27d7ea7..f6d5916 100644
--- a/R/prep.data.AJ.R
+++ b/R/prep.data.AJ.R
@@ -1,163 +1,163 @@
 prep.data.AJ <-
 function(data, states, tr.states){
-    ttime<-c(data$time1,data$Stime)
-    times<-sort(unique(ttime)) # T_ik* indicates the instants of the k transitions
-    
-    
-    # matrix with posible transitions:
-    mat_w<-matrix(FALSE,nrow=3,ncol=3)
-    mat_w[1,2:3]<-TRUE
-    mat_w[2,3]<-TRUE
-    
-    # matrix with all transitions (including censoring transitions):
-    mat_c<-matrix(TRUE,nrow=3,ncol=3)
-    mat_c[2,1]<-FALSE
-    mat_c[3,1:3]<-FALSE
-    
-    colnames(mat_w)<-rownames(mat_w)<-states
-    colnames(mat_c)<-rownames(mat_c)<-states
-    
-    # output states whit contempling censoring (tr 11, tr 22)
-    output.states.c<-lapply(1:dim(mat_c)[2],function(i){
-      rownames(mat_c)[mat_c[i,]==TRUE]})
-    
-    # into states whitout censoring (tr 11, tr 22)
-    into.states<-lapply(1:dim(mat_w)[2],function(i){
-      colnames(mat_w)[mat_w[,i]==TRUE]
-    })
-    
-    
-    # possible transitions (including censoring)
-    
-    to<-c(output.states.c[[1]],output.states.c[[2]])
-    from<-c(rep(tr.states[[1]],length(output.states.c[[1]])),rep(tr.states[[2]],length(output.states.c[[2]])))
-    transitions<-paste("tr",from,to)
-    
-    
-    # risk sets for each non-absorbing state (names)
-    ys <- paste("Y", tr.states) 
-    
-    n=length(data$time1)
-    m=length(times)
-    
-    # number of patientes with time1 = time k-transition (if event1=0 & event=0 --> 1 cens)
-    g11<-integer(m)
-    for(i in 1:m){
-      ss11<-subset(data,data$event1==0 & data$event==0)
-      g11[i]=sum(ss11$time1==times[i])
-    }
-    
-    # number of patients with time1<time k-transition (for transition 1->3)
-    g13<-integer(m)
-    for (i in 1:m){
-      ss13<-subset(data,data$event1==1 & data$event==1 & data$time1==data$Stime)
-      g13[i]=sum(ss13$time1==times[i])
-    }
-    
-    # number of patients with time1<time k-transition (for transition 1->2)
-    g12<-integer(m)
-    for (i in 1:m){
-      ss12<-subset(data,data$time1<data$Stime)
-      g12[i]=sum(ss12$time1==times[i])
-    }
-    
-    # number of patients with time1<time k-transition (for transition 2->2)
-    g22<-integer(m)
-    for (i in 1:m){
-      ss22<-subset(data,data$event1==1 & data$event==0 & data$time1<data$Stime)
-      g22[i]=sum(ss22$Stime==times[i])
-    }
-    
-    # number of patients with time1<time k-transition (for transition 2->3)
-    g23<-integer(m)
-    for (i in 1:m){
-      ss23<-subset(data,data$event1==1 & data$event==1 & data$time1<data$Stime)
-      g23[i]=sum(ss23$Stime==times[i])
-    }
-    
-    ##  matrix of number of transitions:
-    dNs<-cbind(g11,g12,g13,g22,g23)
-    colnames(dNs)<-transitions
-    rownames(dNs)<-times
-    
-    ## matrix of total number of transitions from each non-absorbing state:
-    sum_n1<-rowSums(dNs[,2:3])
-    sum_n2<-dNs[,5]
-    sum_dNs<-cbind(sum_n1,sum_n2)
-    
-    
-    
-    data$start.time<-0
-    
-    initial.state<-integer(length(data[,1]))
-    for(i in 1:length(data[,1])){
-      if(data$start.time[i]==0) initial.state[i]=1
-      else initial.state[i]=2
-    }
-    initial.state<-factor(initial.state,levels=tr.states,labels=tr.states)
-    
-    initial_risk_set<-table(initial.state)
-    
-    
-    
-    
-    
-    i=1
-    n<-initial_risk_set[i]
-    name<-paste("y",i)
-    
-    if(length(into.states[[i]])>0)
-      in.st<-paste("tr",into.states[[i]],i) else in.st<-NULL
-    
-    if(length(output.states.c[[i]])>0)
-      out.st<-paste("tr",i,output.states.c[[i]]) else out.st<-NULL
-    
-    Y1<-c(n,n+cumsum(rowSums(dNs[,in.st,drop=FALSE]))-cumsum(rowSums(dNs[,out.st,drop=FALSE])))
-    
-    i=2
-    n<-initial_risk_set[i]
-    name<-paste("y",i)
-    
-    if(length(into.states[[i]])>0)
-      in.st<-paste("tr",into.states[[i]],i) else in.st<-NULL
-    
-    if(length(output.states.c[[i]])>0)
-      out.st<-paste("tr",i,output.states.c[[i]]) else out.st<-NULL
-    
-    Y2<-c(n,n+cumsum(rowSums(dNs[,in.st,drop=FALSE]))-cumsum(rowSums(dNs[,out.st,drop=FALSE])))
-    
-    Ys<-cbind(Y1[1:m],Y2[1:m])
-    
-    
-    rownames(dNs) <- rownames(sum_dNs) <- rownames(Ys) <- times
-    colnames(dNs) <- transitions
-    colnames(Ys) <- ys
-    colnames(sum_dNs) <- paste("from", tr.states)
-    
-    split_dNs <- strsplit(colnames(dNs), "  ") ## string splits names
-    split_Ys <- strsplit(colnames(Ys), " ")  ## string split names of Ys
-    split_sum_dNs <- strsplit(colnames(sum_dNs), " ") ## string splits names of dNs
-    
-    ## censored columns in dNs, Ys and sum_dNs counting matrix, needed for D-S est
-    cens.dNs.id <- which(sapply(split_dNs, function(x) x[1]=="tr 1 1"|x[1]=="tr 2 2"))
-    cens.Ys.id <- which(sapply(split_Ys, function(x) x[2]%in%tr.states))
-    cens.sum_dNs.id <- which(sapply(split_sum_dNs, function(x) x[2]%in%tr.states))
-    
-    
-    K <- vector(length=nrow(dNs))
-    dN.cens <- rowSums(dNs[, cens.dNs.id, drop=FALSE])
-    Y.cens <- rowSums(Ys[, cens.Ys.id, drop=FALSE]) ## those at risk of being censored
-    N.Y.cens <- ifelse(dN.cens/Y.cens=="NaN", 0, dN.cens/Y.cens)
-    colnames(N.Y.cens) <- NULL
-    H.t <- cumsum(N.Y.cens) ## calculating the hazard
-    k <- exp(-H.t)
-    K <- c(1, k[-length(k)])
-    
-    # weighted for right censoring
-    dNs.w <- dNs/K  ## D-S dNs
-    Ys.w <- Ys/K  ## D-S Ys
-    sum_dNs.w <- sum_dNs/K
-    
-    res <- list(dNs=dNs.w, Ys=Ys.w, sum_dNs=sum_dNs.w, transitions=transitions)
-    return(res)
+  ttime<-c(data$time1,data$Stime)
+  times<-sort(unique(ttime)) # T_ik* indicates the instants of the k transitions
+  
+  
+  # matrix with posible transitions:
+  mat_w<-matrix(FALSE,nrow=3,ncol=3)
+  mat_w[1,2:3]<-TRUE
+  mat_w[2,3]<-TRUE
+  
+  # matrix with all transitions (including censoring transitions):
+  mat_c<-matrix(TRUE,nrow=3,ncol=3)
+  mat_c[2,1]<-FALSE
+  mat_c[3,1:3]<-FALSE
+  
+  colnames(mat_w)<-rownames(mat_w)<-states
+  colnames(mat_c)<-rownames(mat_c)<-states
+  
+  # output states whit contempling censoring (tr 11, tr 22)
+  output.states.c<-lapply(1:dim(mat_c)[2],function(i){
+    rownames(mat_c)[mat_c[i,]==TRUE]})
+  
+  # into states whitout censoring (tr 11, tr 22)
+  into.states<-lapply(1:dim(mat_w)[2],function(i){
+    colnames(mat_w)[mat_w[,i]==TRUE]
+  })
+  
+  
+  # possible transitions (including censoring)
+  
+  to<-c(output.states.c[[1]],output.states.c[[2]])
+  from<-c(rep(tr.states[[1]],length(output.states.c[[1]])),rep(tr.states[[2]],length(output.states.c[[2]])))
+  transitions<-paste("tr",from,to)
+  
+  
+  # risk sets for each non-absorbing state (names)
+  ys <- paste("Y", tr.states) 
+  
+  n=length(data$time1)
+  m=length(times)
+  
+  # number of patientes with time1 = time k-transition (if event1=0 & event=0 --> 1 cens)
+  g11<-integer(m)
+  for(i in 1:m){
+    ss11<-subset(data,data$event1==0 & data$event==0)
+    g11[i]=sum(ss11$time1==times[i])
   }
+  
+  # number of patients with time1<time k-transition (for transition 1->3)
+  g13<-integer(m)
+  for (i in 1:m){
+    ss13<-subset(data,data$event1==1 & data$event==1 & data$time1==data$Stime)
+    g13[i]=sum(ss13$time1==times[i])
+  }
+  
+  # number of patients with time1<time k-transition (for transition 1->2)
+  g12<-integer(m)
+  for (i in 1:m){
+    ss12<-subset(data,data$time1<data$Stime)
+    g12[i]=sum(ss12$time1==times[i])
+  }
+  
+  # number of patients with time1<time k-transition (for transition 2->2)
+  g22<-integer(m)
+  for (i in 1:m){
+    ss22<-subset(data,data$event1==1 & data$event==0 & data$time1<data$Stime)
+    g22[i]=sum(ss22$Stime==times[i])
+  }
+  
+  # number of patients with time1<time k-transition (for transition 2->3)
+  g23<-integer(m)
+  for (i in 1:m){
+    ss23<-subset(data,data$event1==1 & data$event==1 & data$time1<data$Stime)
+    g23[i]=sum(ss23$Stime==times[i])
+  }
+  
+  ##  matrix of number of transitions:
+  dNs<-cbind(g11,g12,g13,g22,g23)
+  colnames(dNs)<-transitions
+  rownames(dNs)<-times
+  
+  ## matrix of total number of transitions from each non-absorbing state:
+  sum_n1<-rowSums(dNs[,2:3])
+  sum_n2<-dNs[,5]
+  sum_dNs<-cbind(sum_n1,sum_n2)
+  
+  
+  
+  data$start.time<-0
+  
+  initial.state<-integer(length(data[,1]))
+  for(i in 1:length(data[,1])){
+    if(data$start.time[i]==0) initial.state[i]=1
+    else initial.state[i]=2
+  }
+  initial.state<-factor(initial.state,levels=tr.states,labels=tr.states)
+  
+  initial_risk_set<-table(initial.state)
+  
+  
+  
+  
+  
+  i=1
+  n<-initial_risk_set[i]
+  name<-paste("y",i)
+  
+  if(length(into.states[[i]])>0)
+    in.st<-paste("tr",into.states[[i]],i) else in.st<-NULL
+  
+  if(length(output.states.c[[i]])>0)
+    out.st<-paste("tr",i,output.states.c[[i]]) else out.st<-NULL
+  
+  Y1<-c(n,n+cumsum(rowSums(dNs[,in.st,drop=FALSE]))-cumsum(rowSums(dNs[,out.st,drop=FALSE])))
+  
+  i=2
+  n<-initial_risk_set[i]
+  name<-paste("y",i)
+  
+  if(length(into.states[[i]])>0)
+    in.st<-paste("tr",into.states[[i]],i) else in.st<-NULL
+  
+  if(length(output.states.c[[i]])>0)
+    out.st<-paste("tr",i,output.states.c[[i]]) else out.st<-NULL
+  
+  Y2<-c(n,n+cumsum(rowSums(dNs[,in.st,drop=FALSE]))-cumsum(rowSums(dNs[,out.st,drop=FALSE])))
+  
+  Ys<-cbind(Y1[1:m],Y2[1:m])
+  
+  
+  rownames(dNs) <- rownames(sum_dNs) <- rownames(Ys) <- times
+  colnames(dNs) <- transitions
+  colnames(Ys) <- ys
+  colnames(sum_dNs) <- paste("from", tr.states)
+  
+  split_dNs <- strsplit(colnames(dNs), "  ") ## string splits names
+  split_Ys <- strsplit(colnames(Ys), " ")  ## string split names of Ys
+  split_sum_dNs <- strsplit(colnames(sum_dNs), " ") ## string splits names of dNs
+  
+  ## censored columns in dNs, Ys and sum_dNs counting matrix, needed for D-S est
+  cens.dNs.id <- which(sapply(split_dNs, function(x) x[1]=="tr 1 1"|x[1]=="tr 2 2"))
+  cens.Ys.id <- which(sapply(split_Ys, function(x) x[2]%in%tr.states))
+  cens.sum_dNs.id <- which(sapply(split_sum_dNs, function(x) x[2]%in%tr.states))
+  
+  
+  K <- vector(length=nrow(dNs))
+  dN.cens <- rowSums(dNs[, cens.dNs.id, drop=FALSE])
+  Y.cens <- rowSums(Ys[, cens.Ys.id, drop=FALSE]) ## those at risk of being censored
+  N.Y.cens <- ifelse(dN.cens/Y.cens=="NaN", 0, dN.cens/Y.cens)
+  colnames(N.Y.cens) <- NULL
+  H.t <- cumsum(N.Y.cens) ## calculating the hazard
+  k <- exp(-H.t)
+  K <- c(1, k[-length(k)])
+  
+  # weighted for right censoring
+  dNs.w <- dNs/K  ## D-S dNs
+  Ys.w <- Ys/K  ## D-S Ys
+  sum_dNs.w <- sum_dNs/K
+  
+  res <- list(dNs=dNs.w, Ys=Ys.w, sum_dNs=sum_dNs.w, transitions=transitions)
+  return(res)
+}
diff --git a/R/prep.data.event.AJ.R b/R/prep.data.event.AJ.R
index e929ef3..54b5d00 100644
--- a/R/prep.data.event.AJ.R
+++ b/R/prep.data.event.AJ.R
@@ -1,25 +1,25 @@
 prep.data.event.AJ <-
 function(dNs, Ys, sum_dNs,states, tr.states) {
-    
-    split_dNs <- strsplit(colnames(dNs), "  ") ## string splits names
-    split_Ys <- strsplit(colnames(Ys), " ")  ## string split names of Ys
-    split_sum_dNs <- strsplit(colnames(sum_dNs), " ") ## string splits names of dNs
-    
-    ## reducing dNs & Ys to just event times & noncens states
-    ##  looks at noncensored columns
-    event.dNs.id <- which(sapply(split_dNs, function(x) x[1]=="tr 1 2" | x[1]=="tr 1 3" | x[1]=="tr 2 3"))
-    ## identifies times where transitions occur
-    event.row.id <- which(apply(dNs[, event.dNs.id, drop=FALSE], 1, function(x) any(x>0)))
-    dNs.event <- dNs[event.row.id, event.dNs.id, drop=FALSE] ## reduces dNs  
-    
-    tr.states.event <- names(which(sapply(tr.states, function(x) length(x)>0))) 
-    event.Ys.id <- which(sapply(split_Ys, function(x) x[2]%in%tr.states.event)) 
-    Ys.event <- Ys[event.row.id, event.Ys.id, drop=FALSE] ## reduces Ys
-    
-    event.sum_dNs.id <- which(sapply(split_sum_dNs, function(x) x[2]%in%states))
-    sum_dNs.event <- sum_dNs[event.row.id, event.sum_dNs.id, drop=FALSE]
-    
-    ans <- list(dNs=dNs.event, Ys=Ys.event, sum_dNs=sum_dNs.event)
-    return(ans)
-    
-  }
+  
+  split_dNs <- strsplit(colnames(dNs), "  ") ## string splits names
+  split_Ys <- strsplit(colnames(Ys), " ")  ## string split names of Ys
+  split_sum_dNs <- strsplit(colnames(sum_dNs), " ") ## string splits names of dNs
+  
+  ## reducing dNs & Ys to just event times & noncens states
+  ##  looks at noncensored columns
+  event.dNs.id <- which(sapply(split_dNs, function(x) x[1]=="tr 1 2" | x[1]=="tr 1 3" | x[1]=="tr 2 3"))
+  ## identifies times where transitions occur
+  event.row.id <- which(apply(dNs[, event.dNs.id, drop=FALSE], 1, function(x) any(x>0)))
+  dNs.event <- dNs[event.row.id, event.dNs.id, drop=FALSE] ## reduces dNs  
+  
+  tr.states.event <- names(which(sapply(tr.states, function(x) length(x)>0))) 
+  event.Ys.id <- which(sapply(split_Ys, function(x) x[2]%in%tr.states.event)) 
+  Ys.event <- Ys[event.row.id, event.Ys.id, drop=FALSE] ## reduces Ys
+  
+  event.sum_dNs.id <- which(sapply(split_sum_dNs, function(x) x[2]%in%states))
+  sum_dNs.event <- sum_dNs[event.row.id, event.sum_dNs.id, drop=FALSE]
+  
+  ans <- list(dNs=dNs.event, Ys=Ys.event, sum_dNs=sum_dNs.event)
+  return(ans)
+  
+}
diff --git a/R/print.TPidm.R b/R/print.TPidm.R
index 748ed50..0f2244a 100644
--- a/R/print.TPidm.R
+++ b/R/print.TPidm.R
@@ -7,64 +7,64 @@ function(x, ...){
     cat("\n")
     
     if(is.null(x$cov)){
-      cat("Parameters:","\n")
-      cat("s=", x$s, "\n")
-      cat("t=", x$t, "\n")
-      cat("Method=", x$method, "\n")
-      cat("CI=", x$CI, "\n")
-      cat("CI transformation=", x$ci.transformation, "\n")
-      cat(paste("Possible transitions:\n", sep=""))
-      print(x$p.trans, row.names = FALSE)
-      
-      
-      cat("\n")
-      if(x$CI==TRUE){
-        if(x$s==0){
-          cat("Occupation probabilities at time t:", "\n", "\n")
-          
-          idm.info<-data.frame(transition=rownames(x$probs),
-                               probs=x$probs[,1],
-                               lower=x$probs[,2],
-                               upper=x$probs[,3],
-                               variance=x$probs[,4])
-          
-        }else{
-          cat("Transition probabilities at time t:", "\n", "\n")
-          
-          idm.info<-data.frame(transition=rownames(x$probs),
-                               probs=x$probs[,1],
-                               lower=x$probs[,2],
-                               upper=x$probs[,3],
-                               variance=x$probs[,4])
-          
+    cat("Parameters:","\n")
+    cat("s=", x$s, "\n")
+    cat("t=", x$t, "\n")
+    cat("Method=", x$method, "\n")
+    cat("CI=", x$CI, "\n")
+    cat("CI transformation=", x$ci.transformation, "\n")
+    cat(paste("Possible transitions:\n", sep=""))
+    print(x$p.trans, row.names = FALSE)
+    
+    
+    cat("\n")
+    if(x$CI==TRUE){
+      if(x$s==0){
+        cat("Occupation probabilities at time t:", "\n", "\n")
+        
+        idm.info<-data.frame(transition=rownames(x$probs),
+                             probs=x$probs[,1],
+                             lower=x$probs[,2],
+                             upper=x$probs[,3],
+                             variance=x$probs[,4])
+        
+      }else{
+        cat("Transition probabilities at time t:", "\n", "\n")
+        
+        idm.info<-data.frame(transition=rownames(x$probs),
+                             probs=x$probs[,1],
+                             lower=x$probs[,2],
+                             upper=x$probs[,3],
+                             variance=x$probs[,4])
+        
+      }
+    }else{
+      if(x$s==0){
+        cat("Occupation probabilities at time t:", "\n", "\n")
+        
+        estimations<-c()
+        for(j in 1:length(x$p.trans)){
+          idx<-unlist(strsplit(x$p.trans[j], " "))
+          estimations[j]<-x$all.probs[length(x$times),,j]
         }
+        
+        idm.info<-data.frame(transition=x$p.trans,
+                             probs=estimations)
+        
       }else{
-        if(x$s==0){
-          cat("Occupation probabilities at time t:", "\n", "\n")
-          
-          estimations<-c()
-          for(j in 1:length(x$p.trans)){
-            idx<-unlist(strsplit(x$p.trans[j], " "))
-            estimations[j]<-x$all.probs[length(x$times),,j]
-          }
-          
-          idm.info<-data.frame(transition=x$p.trans,
-                               probs=estimations)
-          
-        }else{
-          cat("Transition probabilities at time t:", "\n", "\n")
-          
-          estimations<-c()
-          for(j in 1:length(x$p.trans)){
-            idx<-unlist(strsplit(x$p.trans[j], " "))
-            estimations[j]<-x$all.probs[length(x$times),,j]
-          }
-          
-          idm.info<-data.frame(transition=x$p.trans,
-                               probs=estimations)
-          
+        cat("Transition probabilities at time t:", "\n", "\n")
+        
+        estimations<-c()
+        for(j in 1:length(x$p.trans)){
+          idx<-unlist(strsplit(x$p.trans[j], " "))
+          estimations[j]<-x$all.probs[length(x$times),,j]
         }
+        
+        idm.info<-data.frame(transition=x$p.trans,
+                             probs=estimations)
+        
       }
+    }
     } else{
       
       cat("Parameters:","\n")
diff --git a/R/test.nm.R b/R/test.nm.R
index d7e6483..a6b5921 100644
--- a/R/test.nm.R
+++ b/R/test.nm.R
@@ -31,39 +31,39 @@ function(data, s, t="last"){
     
     
     # estimate both methods
-    
-    # method 'NM'
-    # prepare data set and NM method:
-    NM.est<- fun.NM(data, states, s, t, tr.states)
-    
-    # method 'AJ'
-    # start time ==0
-    data$start.time<-0
-    
-    # initial probabilities for each initial state
-    i.state<-integer(length(data[,1]))
-    for(i in 1:length(data[,1])){
-      if(data$start.time[i]==0) i.state[i]=1
-      if(data$start.time[i]==data$time1[i]) i.state[i]=2
-      if(data$start.time[i]==data$Stime[i]) i.state[i]=3
-    }
-    i.state <- factor(i.state,levels=states,labels=states)
-    
-    initial.probs <- prop.table(table(i.state))
-    
-    
-    # prepare data set to compute AJ method:
-    ds.prep.AJ <- prep.data.AJ(data, states, tr.states)
-    
-    # reduces to event times:
-    ds.event.AJ <- prep.data.event.AJ(ds.prep.AJ$dNs, ds.prep.AJ$Ys, ds.prep.AJ$sum_dNs, states, tr.states)
-    event.times <- as.numeric(as.character(rownames(ds.event.AJ$dNs)))
-    
-    # AJ estimator:
-    AJ.est <- fun.AJ(ns,states, ds.event.AJ$dNs, ds.event.AJ$Ys, ds.event.AJ$sum_dNs,
-                     s, t,  event.times, initial.probs)
-    
-    ## plot:
+
+      # method 'NM'
+      # prepare data set and NM method:
+      NM.est<- fun.NM(data, states, s, t, tr.states)
+      
+      # method 'AJ'
+      # start time ==0
+      data$start.time<-0
+      
+      # initial probabilities for each initial state
+      i.state<-integer(length(data[,1]))
+      for(i in 1:length(data[,1])){
+        if(data$start.time[i]==0) i.state[i]=1
+        if(data$start.time[i]==data$time1[i]) i.state[i]=2
+        if(data$start.time[i]==data$Stime[i]) i.state[i]=3
+      }
+      i.state <- factor(i.state,levels=states,labels=states)
+      
+      initial.probs <- prop.table(table(i.state))
+      
+      
+      # prepare data set to compute AJ method:
+      ds.prep.AJ <- prep.data.AJ(data, states, tr.states)
+      
+      # reduces to event times:
+      ds.event.AJ <- prep.data.event.AJ(ds.prep.AJ$dNs, ds.prep.AJ$Ys, ds.prep.AJ$sum_dNs, states, tr.states)
+      event.times <- as.numeric(as.character(rownames(ds.event.AJ$dNs)))
+      
+      # AJ estimator:
+      AJ.est <- fun.AJ(ns,states, ds.event.AJ$dNs, ds.event.AJ$Ys, ds.event.AJ$sum_dNs,
+                       s, t,  event.times, initial.probs)
+      
+      ## plot:
     if(s>0){
       chosen.tr<- c("1 2", "1 3", "2 2")      
       
@@ -86,13 +86,13 @@ function(data, s, t="last"){
       title("Diagnostic plot", outer=TRUE,  line=-1)      
     }
     
-    # results:    
+      # results:    
     cat("Parameters:","\n")
     cat("s=", s, "\n")
     cat("t=", AJ.est$t, "\n")
-    
+        
     cat("\n")
-    if(s==0){
-      cat("Markov assumption is not relevant for the estimation of occupation probabilities (s==0).", "\n", "\n")                
-    }        
+      if(s==0){
+        cat("Markov assumption is not relevant for the estimation of occupation probabilities (s==0).", "\n", "\n")                
+      }        
   }
diff --git a/R/var.AJ.R b/R/var.AJ.R
index f061ee8..a796251 100644
--- a/R/var.AJ.R
+++ b/R/var.AJ.R
@@ -1,84 +1,84 @@
 var.AJ <-
 function(ns, states, dNs.id_tr, Ys.id_tr, sum_dNs.id_tr, TP.AJs, all.I.dA, tr.states){
+  
+  transition.names<-paste(rep(states, ns), rep(states, each=ns))
+  cov.dA<-cov.AJs<-array(0,dim=c(ns^2,ns^2,nrow(dNs.id_tr)),
+                         dimnames=list(rows=transition.names,cols=transition.names, time=rownames(dNs.id_tr)))
+  results.matrix<-array(0,dim=c(ns^2,ns^2)) ## matrix to keep the results to build a cov.AJs matrix
+  colnames(results.matrix)<-rownames(results.matrix)<-transition.names
+  
+  # identity matrices for Kronecker products
+  bl.Id<-diag(ns^2)
+  Id<-diag(ns)
+  
+  vp<-matrix(0,ns,ns)
+  
+  for (i in 1:nrow(dNs.id_tr)) { ## loop through times
     
-    transition.names<-paste(rep(states, ns), rep(states, each=ns))
-    cov.dA<-cov.AJs<-array(0,dim=c(ns^2,ns^2,nrow(dNs.id_tr)),
-                           dimnames=list(rows=transition.names,cols=transition.names, time=rownames(dNs.id_tr)))
-    results.matrix<-array(0,dim=c(ns^2,ns^2)) ## matrix to keep the results to build a cov.AJs matrix
-    colnames(results.matrix)<-rownames(results.matrix)<-transition.names
+    ## VARIANCE OF A-J ( TRANS PROB MATRIX P(s, t) )
+    ## Equation 4.4.20 in Anderson 1993 (The Greenwood type-estimator)
     
-    # identity matrices for Kronecker products
-    bl.Id<-diag(ns^2)
-    Id<-diag(ns)
-    
-    vp<-matrix(0,ns,ns)
-    
-    for (i in 1:nrow(dNs.id_tr)) { ## loop through times
+    ## loop on the blocks (g) - only needed for transitional states (non-absorbing states)
+    for (g in tr.states) {
       
-      ## VARIANCE OF A-J ( TRANS PROB MATRIX P(s, t) )
-      ## Equation 4.4.20 in Anderson 1993 (The Greenwood type-estimator)
+      ## find positioning of g in definition of states
+      id_g <- which(states==g)
+      ## covs: written componentwise, the recursion formula to calculate the variance of transitions in each block (g)
+      covs <- matrix(0, nrow=ns, ncol=ns)
+      colnames(covs) <- rownames(covs) <- states
       
-      ## loop on the blocks (g) - only needed for transitional states (non-absorbing states)
-      for (g in tr.states) {
+      ## loop in the blocks
+      for (j in 1:ns) {
         
-        ## find positioning of g in definition of states
-        id_g <- which(states==g)
-        ## covs: written componentwise, the recursion formula to calculate the variance of transitions in each block (g)
-        covs <- matrix(0, nrow=ns, ncol=ns)
-        colnames(covs) <- rownames(covs) <- states
-        
-        ## loop in the blocks
-        for (j in 1:ns) {
+        ## this just fills in upper diagonal matrix
+        ## use symmetry to fill in rest
+        for (r in j:ns) {
           
-          ## this just fills in upper diagonal matrix
-          ## use symmetry to fill in rest
-          for (r in j:ns) {
-            
-            state.j <- states[j]
-            state.r <- states[r]
-            g.Ys <- paste("Y", g)
-            g.sum_dNs <- paste("from", g)
-            
-            if (Ys.id_tr[i, g.Ys]==0) {  ## if Y_g = 0 then covariance = 0
-              covs[j, r] <- 0
-              next
-            }
+          state.j <- states[j]
+          state.r <- states[r]
+          g.Ys <- paste("Y", g)
+          g.sum_dNs <- paste("from", g)
+          
+          if (Ys.id_tr[i, g.Ys]==0) {  ## if Y_g = 0 then covariance = 0
+            covs[j, r] <- 0
+            next
+          }
+          
+          if (state.j == g & state.r == g) {  ## g==j==r 
+            covs[j, r] <- (Ys.id_tr[i, g.Ys] - sum_dNs.id_tr[i, g.sum_dNs])*sum_dNs.id_tr[i, g.sum_dNs] / Ys.id_tr[i, g.Ys]^3
             
-            if (state.j == g & state.r == g) {  ## g==j==r 
-              covs[j, r] <- (Ys.id_tr[i, g.Ys] - sum_dNs.id_tr[i, g.sum_dNs])*sum_dNs.id_tr[i, g.sum_dNs] / Ys.id_tr[i, g.Ys]^3
-              
-            }  else if (state.j == g & state.r != g) {  ## g!=r
-              name <- paste("tr", g, state.r)
-              if (!name%in%colnames(dNs.id_tr)) next
-              covs[j, r] <- -(Ys.id_tr[i, g.Ys] - sum_dNs.id_tr[i, g.sum_dNs])*dNs.id_tr[i, name] / Ys.id_tr[i, g.Ys]^3
-            } else if (state.j != g & state.r == g) {  ## g!=j
-              name <- paste("tr", g, state.j)
-              if (!name%in%colnames(dNs.id_tr)) next
-              covs[j, r] <- -(Ys.id_tr[i, g.Ys] - sum_dNs.id_tr[i, g.sum_dNs])*dNs.id_tr[i, name]/Ys.id_tr[i, g.Ys]^3
-            } else { ## g!=j and g!=r
-              name.r <- paste("tr", g, state.r)
-              name.j <- paste("tr", g, state.j)
-              if (!(name.j%in%colnames(dNs.id_tr) & name.r%in%colnames(dNs.id_tr))) next
-              covs[j, r] <- (ifelse(j==r, 1, 0)*Ys.id_tr[i, g.Ys] - dNs.id_tr[i, name.j])*dNs.id_tr[i, name.r]/Ys.id_tr[i, g.Ys]^3
-            } ## end of if/else statements
-          } ## end of r loop
-        } ## end of j loop
-        
-        covs[lower.tri(covs)] <- t(covs)[lower.tri(covs)]
-        
-        results.matrix[(seq(1, ns*(ns-1)+1, by=ns) + id_g - 1), (seq(1, ns*(ns-1)+1, by=ns) + id_g - 1)] <- covs
-        
-      }## end of g loop
+          }  else if (state.j == g & state.r != g) {  ## g!=r
+            name <- paste("tr", g, state.r)
+            if (!name%in%colnames(dNs.id_tr)) next
+            covs[j, r] <- -(Ys.id_tr[i, g.Ys] - sum_dNs.id_tr[i, g.sum_dNs])*dNs.id_tr[i, name] / Ys.id_tr[i, g.Ys]^3
+          } else if (state.j != g & state.r == g) {  ## g!=j
+            name <- paste("tr", g, state.j)
+            if (!name%in%colnames(dNs.id_tr)) next
+            covs[j, r] <- -(Ys.id_tr[i, g.Ys] - sum_dNs.id_tr[i, g.sum_dNs])*dNs.id_tr[i, name]/Ys.id_tr[i, g.Ys]^3
+          } else { ## g!=j and g!=r
+            name.r <- paste("tr", g, state.r)
+            name.j <- paste("tr", g, state.j)
+            if (!(name.j%in%colnames(dNs.id_tr) & name.r%in%colnames(dNs.id_tr))) next
+            covs[j, r] <- (ifelse(j==r, 1, 0)*Ys.id_tr[i, g.Ys] - dNs.id_tr[i, name.j])*dNs.id_tr[i, name.r]/Ys.id_tr[i, g.Ys]^3
+          } ## end of if/else statements
+        } ## end of r loop
+      } ## end of j loop
       
-      ## array holding var-cov matrix for I+dA matrix at each time (differential of NA estim)
-      cov.dA[, , i] <- results.matrix
+      covs[lower.tri(covs)] <- t(covs)[lower.tri(covs)]
       
-      if (i==1) {
-        cov.AJs[, , i] <- bl.Id%*% cov.dA[, , i] %*% bl.Id
-      } else {
-        cov.AJs[, , i] <- (t(all.I.dA[, , i]) %x% Id) %*% cov.AJs[, , i-1] %*%((all.I.dA[, , i]) %x% Id) +
-          (Id %x% TP.AJs[, , i-1]) %*% cov.dA[, , i]  %*% (Id%x% t(TP.AJs[, , i-1]))
-      }    
-    } 
-    return(list(TP.AJs=TP.AJs, cov.AJs=cov.AJs))
-  }
+      results.matrix[(seq(1, ns*(ns-1)+1, by=ns) + id_g - 1), (seq(1, ns*(ns-1)+1, by=ns) + id_g - 1)] <- covs
+      
+    }## end of g loop
+    
+    ## array holding var-cov matrix for I+dA matrix at each time (differential of NA estim)
+    cov.dA[, , i] <- results.matrix
+    
+    if (i==1) {
+      cov.AJs[, , i] <- bl.Id%*% cov.dA[, , i] %*% bl.Id
+    } else {
+      cov.AJs[, , i] <- (t(all.I.dA[, , i]) %x% Id) %*% cov.AJs[, , i-1] %*%((all.I.dA[, , i]) %x% Id) +
+        (Id %x% TP.AJs[, , i-1]) %*% cov.dA[, , i]  %*% (Id%x% t(TP.AJs[, , i-1]))
+    }    
+  } 
+  return(list(TP.AJs=TP.AJs, cov.AJs=cov.AJs))
+}
diff --git a/R/var.NM.R b/R/var.NM.R
index 053f213..5fa7248 100644
--- a/R/var.NM.R
+++ b/R/var.NM.R
@@ -48,14 +48,14 @@ function(data, ns, states, tr.states, s, t, probs){
     l.id_tr.s <- length(id_tr.s)
     
     e.times<-times[id_tr.s]
-    
+        
     # Titman variances:
     if (s==0){    
       # subset 1
       set_s1<-subset(data,data$s1==1)
       # status 1:
       n_s1<-nrow(set_s1)
-      
+            
       ### 
       s1_times1<-set_s1$time1   
       s1_Stimes<-set_s1$Stime
@@ -261,8 +261,8 @@ function(data, ns, states, tr.states, s, t, probs){
         A<- sum.var.p12.z.n - (1/n_s1)*sum.var.p12.t.A.n
         B<- sum.var.p12.t.n - (1/n_s1)*sum.var.p12.t.B.n
         
-        S.A <- probs[1,1,i]* A
-        S.B <- (1-probs[1,3,i])* B
+                 S.A <- probs[1,1,i]* A
+                S.B <- (1-probs[1,3,i])* B
         
         SZ_ST <- (S.A-S.B)^2
         
@@ -282,10 +282,10 @@ function(data, ns, states, tr.states, s, t, probs){
       } # end for
       
       ### var(p11)
-      var.p11<-sum.var.p11*((probs[1,1,]/n_s1)^2)
+       var.p11<-sum.var.p11*((probs[1,1,]/n_s1)^2)
       
       ### var(p13)
-      var.p13<-sum.var.p13*(((1-probs[1,3,])/n_s1)^2)
+       var.p13<-sum.var.p13*(((1-probs[1,3,])/n_s1)^2)
       
       ### var(p12)      
       var.p12 <- ((1/n_s1)^2)*sum.var.p12
diff --git a/build/partial.rdb b/build/partial.rdb
index 7aae573..60e1671 100644
Binary files a/build/partial.rdb and b/build/partial.rdb differ
diff --git a/data/colonTP.rda b/data/colonTP.rda
index 364b17c..8d52d96 100644
Binary files a/data/colonTP.rda and b/data/colonTP.rda differ
diff --git a/inst/CITATION b/inst/CITATION
new file mode 100644
index 0000000..87b0f25
--- /dev/null
+++ b/inst/CITATION
@@ -0,0 +1,17 @@
+bibentry(bibtype = "Article",
+  title        = "Estimation of Transition Probabilities for the Illness-Death Model: Package {TP.idm}",
+  author       = c(person(given = "Vanesa",
+                          family = "Balboa"),
+                   person(given = "Jacobo",
+                          family = "de U{\\~n}a-{\\'A}lvarez",
+                          email = "jacobo@uvigo.es")),
+  journal      = "Journal of Statistical Software",
+  year         = "2018",
+  volume       = "83",
+  number       = "10",
+  pages        = "1--19",
+  doi          = "10.18637/jss.v083.i10",
+
+  header       = "To cite TP.idm in publications use:"
+)
+
diff --git a/man/TP.idm-package.Rd b/man/TP.idm-package.Rd
index 53744b3..f83768b 100644
--- a/man/TP.idm-package.Rd
+++ b/man/TP.idm-package.Rd
@@ -30,6 +30,8 @@ Datta S. and Satten G. A. (2001). Validity of the Aalen-Johansen estimators of s
 de Uña-Álvarez J. and Meira-Machado L. (2015). Nonparametric estimation of transition probabilities in a non-Markov illness-death model: a comparative study. \emph{Biometrics} \bold{71}, 364--375.
 
 Meira-Machado L. F., de Uña-Álvarez J. and Cadarso-Suárez C. (2006). Nonparametric estimation of transition probabilities in a non-Markov illness-death model. \emph{Lifetime Data Annal} \bold{12(3)}, 325--344.
+
+Balboa-Barreiro V., de Uña-Álvarez J. (2018). Estimation of Transition Probabilities for the Illness-Death Model: Package \code{TP.idm}. \emph{Journal of Statistical Software} \bold{83(10)}, 1--19.
 }
 
 \keyword{ package }
diff --git a/man/TPidm.Rd b/man/TPidm.Rd
index 547f464..0257267 100644
--- a/man/TPidm.Rd
+++ b/man/TPidm.Rd
@@ -40,7 +40,7 @@ ci.transformation = "linear", method = "NM")
 }
 }
 \details{
-  If s = 0 this function TPidm reports the state occupation probabilities at time t. For s > 0, the transition probabilities are provided. The default method \dQuote{NM} computes the estimator proposed by de Uña-Álvarez and Meira-Machado (2015), which is consistent regardless the Markov assumption. To fit a Markovian transition probability matrix, use the \dQuote{AJ} method. Note that the Aalen-Johansen occupation probabilities are consistent even when the process is non-Markov (Datta and Satten, 2001), but this is not true for the Aalen-Johansen transition probabilities (Meira-Machado et al., 2006). If \code{CI} is \code{TRUE} (default), TPidm calculates the estimated variance and the confidence intervals by using the plug-in variance (method = \dQuote{NM}), or the standard formulas for Markovian processes (method = \dQuote{AJ}, see Andersen et al., 1993).
+  If s = 0 this function TPidm reports the state occupation probabilities at time t. For s > 0, the transition probabilities are provided. The default method \dQuote{NM} computes the estimator proposed by de Uña-Álvarez and Meira-Machado (2015), which is consistent regardless the Markov assumption. To fit a Markovian transition probability matrix, use the \dQuote{AJ} method. Note that the Aalen-Johansen occupation probabilities are consistent even when the process is non-Markov (Datta and Satten, 2001), but this is not true for the Aalen-Johansen transition probabilities (Meira-Machado et al., 2006). If \code{CI} is \code{TRUE} (default), TPidm calculates the estimated variance and the confidence intervals by using the plug-in variance described in Balboa and de Uña-Álvarez (2018) (method = \dQuote{NM}), or the standard formulas for Markovian processes (method = \dQuote{AJ}, see Andersen et al., 1993).
   }
 \value{
   An object of class \sQuote{TPidm}. It is a list containing the following objects:
@@ -68,6 +68,8 @@ If \code{cov} is specified, the main function provides estimations of transition
  de Uña-Álvarez J. and Meira-Machado L. (2015). Nonparametric estimation of transition probabilities in a non-Markov illness-death model: a comparative study. \emph{Biometrics} \bold{71}, 364--375.
  
  Meira-Machado L. F., de Uña-Álvarez J. and Cadarso-Suárez C. (2006). Nonparametric estimation of transition probabilities in a non-Markov illness-death model. \emph{Lifetime Data Anal} \bold{12(3)}, 325--344.
+ 
+ Balboa-Barreiro V., de Uña-Álvarez J. (2018). Estimation of Transition Probabilities for the Illness-Death Model: Package \code{TP.idm}. \emph{Journal of Statistical Software} \bold{83(10)}, 1--19.
 }
 %%\author{
 %%  ~~who you are~~
diff --git a/man/colonTP.Rd b/man/colonTP.Rd
index 3439493..6554f2c 100644
--- a/man/colonTP.Rd
+++ b/man/colonTP.Rd
@@ -30,7 +30,7 @@
   }
 }
 \note{
-  The study is originally described in Laurie et al. (1989). The main report is found in Moertel et al. (1990). This data set is closet to that of the final report in Moertel et al. (1991). A version of the data with less follow-up time was used in the paper by Lin (1994). This data set has been used for illustration purposes in a number of papers, including Meira-Machado et al. (2014) and de Uña-Álvarez and Meira-Machado (2015).
+  The study is originally described in Laurie et al. (1989). The main report is found in Moertel et al. (1990). This data set is closet to that of the final report in Moertel et al. (1991). A version of the data with less follow-up time was used in the paper by Lin (1994). This data set has been used for illustration purposes in a number of papers, including Meira-Machado et al. (2015) and de Uña-Álvarez and Meira-Machado (2015).
 }
 %%\source{
 %%  ~~ reference to a publication or URL from which the data were obtained ~~
@@ -42,7 +42,7 @@
   
   Lin D. Y. (1994). Cox regression analysis of multivariate failure time data: the marginal approach. \emph{Statistics in Medicine} \bold{13(21)}, 2233--2247.
   
-  Meira-Machado L., de Uña-Álvarez J. and Datta S. (2014). Nonparametric estimation of conditional transition probabilities in a non-Markov illness-death model. Computational Statistics, in press. <DOI:10.1007/s00180-014-0538-6>.
+  Meira-Machado L., de Uña-Álvarez J. and Datta S. (2015). Nonparametric estimation of conditional transition probabilities in a non-Markov illness-death model. Computational Statistics \bold{30}, 377--397.
   
   Moertel C. G., Fleming T. R., MacDonald J. S., Haller D. G., Laurie J. A., Goodman P. J., Ungerleider J. S., Emerson W. A., Tormey D. C., Glick J. H., Veeder M. H. and Maillard J. A. (1990). Levamisole and fluorouracil for adjuvant therapy of resected colon carcinoma. \emph{New England Journal of Medicine} \bold{322(6)}, 352--358.
   
diff --git a/src/init.c b/src/init.c
new file mode 100644
index 0000000..9fc6acd
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,20 @@
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+/* FIXME: 
+  Check these declarations against the C/Fortran source code.
+*/
+  
+  /* .C calls */
+  extern void WeightsKaplanMeierSort(void *, void *, void *, void *, void *);
+
+static const R_CMethodDef CEntries[] = {
+  {"WeightsKaplanMeierSort", (DL_FUNC) &WeightsKaplanMeierSort, 5},
+  {NULL, NULL, 0}
+};
+
+void R_init_TP_idm(DllInfo *dll)
+{
+  R_registerRoutines(dll, CEntries, NULL, NULL, NULL);
+  R_useDynamicSymbols(dll, FALSE);
+}