depmixS An Package for Hidden Markov Models Ingmar Visser University of Amsterdam Maarten - PDF document

depmixS4:AnRPackageforHiddenMarkovModelsIngmarVisserUniversityofAmsterdamMaartenSpeekenbrinkUniversityCollegeLondon AbstractThisintroductiontotheRpackagedepmixS4isa(slightly)modiedversionofVisserandSpeekenbrink(2010),publishedintheJournalofStatisticalSoftware.PleaserefertothatarticlewhenusingdepmixS4.Thecurrentversionis1.4-2;theversionhistoryandchangescanbefoundintheNEWSleofthepackage.Below,themajorversionsarelistedalongwiththemostnoteworthychanges.depmixS4implementsageneralframeworkfordeningandestimatingdependentmix-turemodelsintheRprogramminglanguage.ThisincludesstandardMarkovmodels,la-tent/hiddenMarkovmodels,andlatentclassandnitemixturedistributionmodels.Themodelscanbettedonmixedmultivariatedatawithdistributionsfromtheglmfamily,the(logistic)multinomial,orthemultivariatenormaldistribution.Otherdistributionscanbeaddedeasily,andanexampleisprovidedwiththeexgausdistribution.Parametersareestimatedbytheexpectation-maximization(EM)algorithmor,when(linear)con-straintsareimposedontheparameters,bydirectnumericaloptimizationwiththeRsolnporRdonlp2routines.Keywords:hiddenMarkovmodel,dependentmixturemodel,mixturemodel,constraints. VersionhistorySeetheNEWSleforcompletelistingsofchanges;hereonlythemajorchangesarementioned.1.3-0Addedoption`classication'likelihoodtoEM;modeloutputisnowpretty-printedandparametersaregivenpropernames;thetfunctionhasgainedargumentsfornecontrolofusingRsolnpandRdonlp2.1.2-0Addedsupportformissingdata,seesection2.3.1.1-0SpeedimprovementsduetowritingthemainloopinCcode.1.0-0Firstreleasewiththisvignette,amodiedversionofthepaperintheJournalofStatisticalSoftware.0.1-0FirstversiononCRAN.1.IntroductionMarkovandlatentMarkovmodelsarefrequentlyusedinthesocialsciences,indierentareasandapplications.Inpsychology,theyareusedformodellinglearningprocesses;seeWickens 2depmixS4:AnRPackageforHiddenMarkovModels(1982),foranoverview,ande.g.,Schmittmann,Visser,andRaijmakers(2006),forarecentapplication.Ineconomics,latentMarkovmodelsareso-calledregimeswitchingmodels(seee.g.,Kim1994andGhysels1994).Furtherapplicationsincludespeechrecognition(Rabiner1989),EEGanalysis(RainerandMiller2000),andgenetics(Krogh1998).Intheselatterareasofapplication,latentMarkovmodelsareusuallyreferredtoashiddenMarkovmodels.SeeforexampleFruhwirth-Schnatter(2006)foranoverviewofhiddenMarkovmodelswithextensions.Furtherexamplesofapplicationscanbefoundine.g.,Cappe,Moulines,andRyden(2005,Chapter1).AmoregentleintroductionintohiddenMarkovmodelswithapplicationsisthebookbyZucchiniandMacDonald(2009).ThedepmixS4packagewasmotivatedbythefactthatwhileMarkovmodelsareusedcom-monlyinthesocialsciences,nocomprehensivepackagewasavailableforttingsuchmodels.ExistingsoftwareforestimatingMarkovianmodelsincludePanmark(vandePol,Langeheine,andJong1996),andforlatentclassmodelsLatentGold(VermuntandMagidson2003).Theseprogramslackanumberofimportantfeatures,besidesnotbeingfreelyavailable.TherearecurrentlysomepackagesinRthathandlehiddenMarkovmodelsbuttheylackanumberoffeaturesthatweneededinourresearch.Inparticular,depmixS4wasdesignedtomeetthefollowinggoals:1.tobeabletoestimateparameterssubjecttogenerallinear(in)equalityconstraints;2.tobeabletottransitionmodelswithcovariates,i.e.,tohavetime-dependenttransitionmatrices;3.tobeabletoincludecovariatesinthepriororinitialstateprobabilities;4.tobeeasilyextensible,inparticular,toallowuserstoeasilyaddnewuni-ormultivariateresponsedistributionsandnewtransitionmodels,e.g.,continuoustimeobservationmodels.AlthoughdepmixS4wasdesignedtodealwithlongitudinalortimeseriesdata,forsayT�100,itcanalsohandlethelimitcasewhenT=1.Inthiscase,therearenotimedependenciesbetweenobserveddataandthemodelreducestoanitemixtureorlatentclassmodel.Whiletherearespecializedpackagestodealwithmixturedata,asfarasweknowthesedonotallowtheinclusionofcovariatesonthepriorprobabilitiesofclassmembership.ThepossibilitytoestimatetheeectsofcovariatesonpriorandtransitionprobabilitiesisadistinguishingfeatureofdepmixS4.InSection2,weprovideanoutlineofthemodelandlikelihoodequations.ThedepmixS4packageisimplementedusingtheRsystemforstatisticalcomputing(RDe-velopmentCoreTeam2010)andisavailablefromtheComprehensiveRArchiveNetworkathttp://CRAN.R-project.org/package=depmixS4.2.ThedependentmixturemodelThedataconsideredherehavethegeneralformO1:T=(O11;:::;Om1,O12;:::;Om2,...,O1T;:::;OmT)foranm-variatetimeseriesoflengthT.Inthefollowing,weuseOtasshort-handforO1t;:::;Omt.Asanexample,consideratimeseriesofresponsesgeneratedbyasingleparticipantinapsychologicalresponsetimeexperiment.Thedataconsistsofthreevariables:responsetime,responseaccuracy,andacovariatewhichisapay-ovariablere ectingthe IngmarVisser,MaartenSpeekenbrink3 Figure1:Responsetimes(rt),accuracy(corr)andpay-ovalues(Pacc)fortherstseriesofresponsesindatasetspeed.relativerewardforspeededand/oraccurateresponding.Thesevariablesaremeasuredon168,134and137occasionsrespectively(therstseriesof168trialsisplottedinFigure1).ThesedataaremorefullydescribedinDutilh,Wagenmakers,Visser,andvanderMaas(2011),andinthenextsectionanumberofexamplemodelsforthesedataisdescribed.ThelatentMarkovmodelisusuallyassociatedwithdataofthistype,inparticularformulti-nomiallydistributedresponses.However,commonlyemployedestimationprocedures(e.g.,vandePoletal.1996),arenotsuitableforlongtimeseriesduetounder owproblems.Incontrast,thehiddenMarkovmodelistypicallyonlyusedfor`long'univariatetimeseries(Cappeetal.2005,Chapter1).Weusetheterm\dependentmixturemodel"becauseoneoftheauthors(IngmarVisser)thoughtitwastimeforanewnametorelatethesemodels1.Thefundamentalassumptionofadependentmixturemodelisthatatanytimepoint,theobservationsaredistributedasamixturewithncomponents(orstates),andthattime-dependenciesbetweentheobservationsareduetotime-dependenciesbetweenthemixturecomponents(i.e.,transitionprobabilitiesbetweenthecomponents).Theselatterdependenciesareassumedtofollowarst-orderMarkovprocess.Inthemodelsweareconsideringhere,themixturedistributions,theinitialmixtureprobabilitiesandtransitionprobabilitiescanalldependoncovariateszt.Inadependentmixturemodel,thejointlikelihoodofobservationsO1:TandlatentstatesS1:T=(S1;:::;ST),givenmodelparametersandcovariatesz1:T=(z1;:::;zT),canbe 1OnlylaterhefoundoutthatLerouxandPuterman(1992)alreadycoinedthetermdependentmixturemodelsinanapplicationwithhiddenMarkovmixturesofPoissoncountdata. 4depmixS4:AnRPackageforHiddenMarkovModelswrittenas:P(O1:T;S1:Tj;z1:T)=i(z1)bS1(O1jz1)T�1Yt=1aij(zt)bSt(Ot+1jzt+1);(1)wherewehavethefollowingelements:1.StisanelementofS=f1:::ng,asetofnlatentclassesorstates.2.i(z1)=P(S1=ijz1),givingtheprobabilityofclass/stateiattimet=1withcovariatez1.3.aij(zt)=P(St+1=jjSt=i;zt),providestheprobabilityofatransitionfromstateitostatejwithcovariatezt.4.bStisavectorofobservationdensitiesbkj(zt)=P(OktjSt=j;zt)thatprovidetheconditionaldensitiesofobservationsOktassociatedwithlatentclass/statejandcovariatezt,j=1;:::;n,k=1;:::;m.Fortheexampledataabove,bkjcouldbeaGaussiandistributionfunctionfortheresponsetimevariable,andaBernoullidistributionfortheaccuracyvariable.Inthemodelsconsideredhere,boththetransitionprobabilityfunctionsaijandtheinitialstateprobabilityfunctionsmaydependoncovariatesaswellastheresponsedistributionsbkj.2.1.LikelihoodToobtainmaximumlikelihoodestimatesofthemodelparameters,weneedthemarginallikelihoodoftheobservations.ForhiddenMarkovmodels,thismarginal(log-)likelihoodisusuallycomputedbytheso-calledforward-backwardalgorithm(BaumandPetrie1966;Rabiner1989),orratherbytheforwardpartofthisalgorithm.LystigandHughes(2002)changedtheforwardalgorithminsuchawayastoallowcomputingthegradientsofthelog-likelihoodatthesametime.Theystartbyrewritingthelikelihoodasfollows(foreaseofexpositionthedependenceonthemodelparametersandcovariatesisdroppedhere):LT=P(O1:T)=TYt=1P(OtjO1:(t�1));(2)whereP(O1jO0):=P(O1).Notethatforasimple,i.e.,observed,Markovchaintheseproba-bilitiesreducetoP(OtjO1;:::;Ot�1)=P(OtjOt�1).Thelog-likelihoodcannowbeexpressedas:lT=TXt=1log[P(OtjO1:(t�1)]:(3)Tocomputethelog-likelihood,LystigandHughes(2002)denethefollowing(forward)re-cursion:1(j):=P(O1;S1=j)=jbj(O1)(4)t(j):=P(Ot;St=jjO1:(t�1))=nXi=1[t�1(i)aijbj(Ot)](t�1)�1;(5) IngmarVisser,MaartenSpeekenbrink5wheret=Pni=1t(i).Combiningt=P(OtjO1:(t�1)),andequation(3)givesthefollowingexpressionforthelog-likelihood:lT=TXt=1logt:(6)2.2.ParameterestimationParametersareestimatedindepmixS4usingtheexpectation-maximization(EM)algorithmorthroughtheuseofageneralNewton-Raphsonoptimizer.IntheEMalgorithm,parametersareestimatedbyiterativelymaximizingtheexpectedjointlog-likelihoodoftheparametersgiventheobservationsandstates.Let=(1;2;3)bethegeneralparametervectorconsistingofthreesubvectorswithparametersforthepriormodel,transitionmodel,andresponsemodelsrespectively.Thejointlog-likelihoodcanbewrittenas:logP(O1:T;S1:Tjz1:T;)=logP(S1jz1;1)+TXt=2logP(StjSt�1;zt�1;2)+TXt=1logP(OtjSt;zt;3)(7)ThislikelihooddependsontheunobservedstatesS1:T.IntheExpectationstep,wereplacethesewiththeirexpectedvaluesgivenasetof(initial)parameters0=(01;02;03)andobservationsO1:T.Theexpectedlog-likelihood:Q(;0)=E0(logP(O1:T;S1:TjO1:T;z1:T;));(8)canbewrittenas:Q(;0)=nXj=1 1(j)logP(S1=jjz1;1)+TXt=2nXj=1nXk=1t(j;k)logP(St=kjSt�1=j;zt�1;2)+TXt=1nXj=1mXk=1 t(j)logP(OktjSt=j;zt;3);(9)wheretheexpectedvaluest(j;k)=P(St=k;St�1=jjO1:T;z1:T;0)and t(j)=P(St=jjO1:T;z1:T;0)canbecomputedeectivelybytheforward-backwardalgorithm(seee.g.,Rabiner1989).TheMaximizationstepconsistsofthemaximizationof(9)for.Astherighthandsideof(9)consistsofthreeseparateparts,wecanmaximizeseparatelyfor1,2and3.Incommonmodels,maximizationfor1and2isperformedbythennet.defaultroutineinthennetpackage(VenablesandRipley2002),andmaximizationfor3bythestandardglmroutine.Notethatforthelattermaximization,theexpectedvalues t(j)areusedaspriorweightsoftheobservationsOkt.TheEMalgorithmhoweverhassomedrawbacks.First,itcanbeslowtoconvergetowardstheendofoptimization.Second,applyingconstraintstoparameterscanbeproblematic;in 6depmixS4:AnRPackageforHiddenMarkovModelsparticular,EMcanleadtowrongparameterestimateswhenapplyingconstraints.Hence,indepmixS4,EMisusedbydefaultinunconstrainedmodels,butotherwise,directoptimizationisused.TwooptionsareavailablefordirectoptimizationusingpackageRsolnp(GhalanosandTheul2010;Ye1987),orRdonlp2(Tamura2009;Spellucci2002).Bothpackagescanhandlegenerallinear(in)equalityconstraints(andoptionallyalsonon-linearconstraints).2.3.MissingdataMissingdatacanbedealtwithstraightforwardlyincomputingthelikelihoodusingtheforwardrecursioninEquations(4{5).AssumewehaveobservedO1:(t�1)butthatobservationOtismissing.Thekeyideathat,inthiscase,thelteringdistribution,theprobabilitiest,shouldbeidenticaltothestatepredictiondistribution,asthereisnoadditionalinformationtoestimatethecurrentstate.Thus,theforwardvariablestarenowcomputedas:t(i)=P(St=ijO1:(t�1))(10)=nXj=1t�1(j)P(St=ijSt�1=j):(11)Forlaterobservations,wecanthenusethislatterequationagain,realizingthatthelteringdistributionistechnicallye.g.P(St+1jO1:(t�1);t+1).Computationally,theeasiestwaytoimplementthisistosimplysetb(OtjSt)=1ifOtismissing.3.UsingdepmixS4TwostepsareinvolvedinusingdepmixS4whichareillustratedbelowwithexamples:1.modelspecicationwithfunctiondepmix(orwithmixforlatentclassandnitemixturemodels,seeexamplebelowonaddingcovariatestopriorprobabilities);2.modelttingwithfunctionfit.Wehaveseparatedthestagesofmodelspecicationandmodelttingbecausettinglargemodelscanbefairlytime-consuminganditishenceusefultobeabletocheckthemodelspecicationbeforeactuallyttingthemodel.3.1.Exampledata:speedThroughoutthisarticleadatasetcalledspeedisused.Asalreadyindicatedintheintroduc-tion,itconsistsofthreetimeserieswiththreevariables:responsetimert,accuracycorr,andacovariate,Pacc,whichdenestherelativepay-oforspeededversusaccurateresponding.Beforedescribingsomeofthemodelsthatarettedtothesedata,weprovideabriefsketchofthereasonsforgatheringthesedataintherstplace.Responsetimesareaverycommondependentvariableinpsychologicalexperimentsandhenceformthebasisforinferenceaboutmanypsychologicalprocesses.Apotentialthreattosuchinferencebasedonresponsetimesisformedbythespeed-accuracytrade-o:dierentparticipantsinanexperimentmayresponddierentlytotypicalinstructionsto`respondasfastandaccurateaspossible'.Apopularmodelwhichtakesthespeed-accuracytrade-o IngmarVisser,MaartenSpeekenbrink7intoaccountisthediusionmodel(Ratcli1978),whichhasproventoprovideaccuratedescriptionsofresponsetimesinmanydierentsettings.Onepotentialproblemwiththediusionmodelisthatitpredictsacontinuoustrade-obetweenspeedandaccuracyofresponding,i.e.,whenparticipantsarepressedtorespondfasterandfaster,thediusionmodelpredictsthatthiswouldleadtoagradualdecreaseinaccuracy.Thespeeddatasetthatweanalyzebelowwasgatheredtotestthishypothesisversusthealternativehypothesisstatingthatthereisasuddentransitionfromslowandaccuraterespondingtofastrespondingatchancelevel.Ateachtrialoftheexperiment,theparticipantisshownthecurrentsettingoftherelativerewardforspeedversusaccuracy.ThebottompanelofFigure1showsthevaluesofthisvariable.Theexperimentwasdesignedtoinvestigatewhatwouldhappenwhenthisrewardvariablechangesfromrewardforaccuracyonlytorewardforspeedonly.ThespeeddatathatweanalyseherearefromparticipantAinExperiment1inDutilhetal.(2011),whoprovideacompletedescriptionoftheexperimentandtherelevanttheoreticalbackground.Thecentralquestionregardingthisdataiswhetheritisindeedbestdescribedbytwomodesofrespondingratherthanasinglemodeofrespondingwithacontinuoustrade-obetweenspeedandaccuracy.Thehallmarkofadiscontinuitybetweenslowversusspeededrespondingisthatswitchingbetweenthetwomodesisasymmetric(seee.g.VanderMaasandMolenaar1992,foratheoreticalunderpinningofthisclaim).ThefithelppageofdepmixS4providesanumberofexamplesinwhichtheasymmetryoftheswitchingprocessistested;thoseexamplesandothercandidatemodelsarediscussedatlengthinVisser,Raijmakers,andVanderMaas(2009).3.2.AsimplemodelAdependentmixturemodelisdenedbythenumberofstatesandtheinitialstate,statetransition,andresponsedistributionfunctions.Adependentmixturemodelcanbecreatedwiththedepmixfunctionasfollows:R�library("depmixS4")R�data("speed")R�set.seed(1)R�mod-depmix(response=rt~1,data=speed,nstates=2,+trstart=runif(4))TherstlineofcodeloadsthedepmixS4packageanddata(speed)loadsthespeeddataset.Thelineset.seed(1)isnecessarytogetstartingvaluesthatwillresultintherightmodel,seemoreonstartingvaluesbelow.Thecalltodepmixspeciesthemodelwithanumberofarguments.Theresponseargumentisusedtospecifytheresponsevariable,possiblywithcovariates,inthefamiliarformatusingformulaesuchasinlmorglmmodels.Thesecondargument,data=speed,providesthedata.frameinwhichthevariablesfromresponseareinterpreted.Third,thenumberofstatesofthemodelisgivenbynstates=2.Startingvalues.Notealsothatstartvaluesforthetransitionparametersareprovidedinthiscallbythetrstartargument.Startingvaluesforparameterscanbeprovidedusing 8depmixS4:AnRPackageforHiddenMarkovModelsthreearguments:instartfortheparametersrelatingtothepriorprobabilities,trstart,relatingthetransitionmodels,andrespstartfortheparametersoftheresponsemodels.Notethatthestartingvaluesfortheinitialandtransitionmodelsaswellasmultinomiallogitresponsemodelsareinterpretedasprobabilities,andinternallyconvertedtomultinomiallogitparameters(ifalogitlinkfunctionisused).StartvaluescanalsobegeneratedrandomlywithintheEMalgorithmbygeneratingrandomuniformvaluesforthevaluesof tin(9)atiteration0.Automaticgenerationofstartingvaluesisnotavailableforconstrainedmodels(seebelow).Inthecalltofitbelow,theargumentemc=em.control(rand=FALSE)controlstheEMalgorithmandspeciesthatrandomstartvaluesshouldnotbegenerated2.Fittingthemodelandprintingresults.Thedepmixfunctionreturnsanobjectofclass`depmix'whichcontainsthemodelspecication,andnotattedmodel.Hence,themodelneedstobettedbycallingfit:R�fm-fit(mod,emc=em.control(rand=FALSE))convergedatiteration68withlogLik:-88.73Thefitfunctionreturnsanobjectofclass`depmix.fitted'whichextendsthe`depmix'class,addingconvergenceinformation(andinformationaboutconstraintsifthesewereapplied,seebelow).Theclasshasprintandsummarymethodstoseetheresults.Theprintmethodprovidesinformationonconvergence,thelog-likelihoodandtheAICandBICvalues:RကfmConvergenceinfo:Loglikelihoodconvergedtowithintol.(relativechange)'logLik.'-88.73(df=7)AIC:191.5BIC:220.1ThesestatisticscanalsobeextractedusinglogLik,AICandBIC,respectively.Bycompari-son,a1-statemodelforthesedata,i.e.,assumingthereisnomixture,hasalog-likelihoodof�305:33,and614.66,and622.83fortheAICandBICrespectively.Hence,the2-statemodeltsthedatamuchbetterthanthe1-statemodel.Notethatthe1-statemodelcanbespec-iedusingmod-depmix(rt~1,data=speed,nstates=1),althoughthismodelistrivialasitwillsimplyreturnthemeanandstandarddeviationofthertvariable.Thesummarymethodoffittedmodelsprovidestheparameterestimates,rstforthepriorprobabilitiesmodel,secondforthetransitionmodels,andthirdfortheresponsemodels.Rကsummary(fm)Initialstateprobabilitiesmodelpr1pr210 2Asofversion1.0-1,thedefaultishaverandomparameterinitializationwhenusingtheEMalgorithm. IngmarVisser,MaartenSpeekenbrink9TransitionmatrixtoS1toS2fromS10.9160.084fromS20.1160.884ResponseparametersResp1:gaussianRe1.(Intercept)Re1.sdSt16.3850.244St25.5100.192Sincenofurtherargumentswerespecied,theinitialstate,statetransitionandresponsedistributionsweresettotheirdefaults(multinomialdistributionsforthersttwo,andaGaussiandistributionfortheresponse).Theresultingmodelindicatestwowell-separatedstates,onewithslowandthesecondwithfastresponses.Thetransitionprobabilitiesindicateratherstablestates,i.e.,theprobabilityofremainingineitherofthestatesisaround0.9.Theinitialstateprobabilityestimatesindicatethatstate1isthestartingstatefortheprocess,withanegligibleprobabilityofstartinginstate2.3.3.CovariatesontransitionparametersBydefault,thetransitionprobabilitiesandtheinitialstateprobabilitiesareparameterizedusingamultinomialmodelwithanidentitylinkfunction.Usingamultinomiallogisticmodelallowsonetoincludecovariatesontheinitialstateandtransitionprobabilities.Inthiscase,eachrowofthetransitionmatrixisparameterizedbyabaselinecategorylogisticmultinomial,meaningthattheparameterforthebasecategoryisxedatzero(seeAgresti2002,p.267.,formultinomiallogisticmodelsandvariousparameterizations).Thedefaultbaselinecate-goryistherststate.Hence,forexample,fora3-statemodel,theinitialstateprobabilitymodelwouldhavethreeparametersofwhichtherstisxedatzeroandtheothertwoarefreelyestimated.Chung,Walls,andPark(2007)discussarelatedlatenttransitionmodelforrepeatedmeasurementdata(T=2)usinglogisticregressiononthetransitionparameters;theyrelyonBayesianmethodsofestimation.Covariatesonthetransitionprobabilitiescanbespeciedusingaone-sidedformulaasinthefollowingexample:R�set.seed(1)R�mod-depmix(rt~1,data=speed,nstates=2,family=gaussian(),+transition=~scale(Pacc),instart=runif(2))Rကfm-fit(mod,verbose=FALSE,emc=em.control(rand=FALSE))convergedatiteration42withlogLik:-44.2Notetheuseofverbose=FALSEtosuppressprintingofinformationontheiterationsofthettingprocess.Applyingsummarytothefittedobjectgives(onlytransitionmodelsprintedherebyusingargumentwhich):Rကsummary(fm,which="transition") 10depmixS4:AnRPackageforHiddenMarkovModelsTransitionmodelforstate(component)1Modeloftypemultinomial(mlogit),formula:~scale(Pacc)Coefficients:St1St2(Intercept)0-0.9518scale(Pacc)01.3924Probalitiesatzerovaluesofthecovariates.0.72150.2785Transitionmodelforstate(component)2Modeloftypemultinomial(mlogit),formula:~scale(Pacc)Coefficients:St1St2(Intercept)02.472scale(Pacc)03.581Probalitiesatzerovaluesofthecovariates.0.077880.9221Thesummaryprovidesallparametersofthemodel,alsothe(redundant)zeroesforthebase-linecategoryinthemultinomialmodel.Thesummaryalsoprintsthetransitionprobabilitiesatthezerovalueofthecovariate.Notethatscalingofthecovariateisusefulinthisregardasitmakesinterpretationoftheseinterceptprobabilitieseasier.3.4.MultivariatedataMultivariatedatacanbemodelledbyprovidingalistofformulaeaswellasalistoffamilyobjectsforthedistributionsofthevariousresponses.InaboveexampleswehaveonlyusedtheresponsetimeswhichweremodelledasaGaussiandistribution.Theaccuracyvariableinthespeeddatacanbemodelledwithamultinomialbyspecifyingthefollowing:R�set.seed(1)R�mod-depmix(list(rt~1,corr~1),data=speed,nstates=2,+family=list(gaussian(),multinomial("identity")),+transition=~scale(Pacc),instart=runif(2))Rကfm-fit(mod,verbose=FALSE,emc=em.control(rand=FALSE))convergedatiteration31withlogLik:-255.5Thisprovidesthefollowingttedmodelparameters(onlytheresponseparametersaregivenhere):Rကsummary(fm,which="response")ResponseparametersResp1:gaussianResp2:multinomialRe1.(Intercept)Re1.sdRe2.incRe2.corSt15.5170.1970.4750.525St26.3910.2390.0980.902 IngmarVisser,MaartenSpeekenbrink11Ascanbeseen,state1hasfastresponsetimesandaccuracyisapproximatelyatchancelevel(.4747),whereasstate2correspondswithslowerrespondingathigheraccuracylevels(.9021).Notethatbyspecifyingmultivariateobservationsintermsofalist,thevariablesarecon-sideredconditionallyindependent(giventhestates).Conditionallydependentvariablesmustbehandledasasingleelementinthelist.Eectively,thismeansspecifyingamultivari-ateresponsemodel.Currently,depmixS4hasonemultivariateresponsemodelwhichisformultivariatenormalvariables.3.5.FixingandconstrainingparametersUsingpackageRsolnp(GhalanosandTheul2010)orRdonlp2(Tamura2009),parametersmaybettedsubjecttogenerallinear(in-)equalityconstraints.Constrainingandxingparametersisdoneusingtheconpatargumenttothefitfunction,whichspeciesforeachparameterinthemodelwhetherit'sxed(0)orfree(1orhigher).Equalityconstraintscanbeimposedbygivingtwoparametersthesamenumberintheconpatvector.Whenonlyxedvaluesarerequired,thefixedargumentcanbeusedinsteadofconpat,withzeroesforxedparametersandothervalues(e.g.,ones)fornon-xedparameters.Fittingthemodelssubjecttotheseconstraintsishandledbytheoptimizationroutinesolnpor,optionally,bydonlp2.Tobeabletoconstructtheconpatand/orfixedvectorsoneneedsthecorrectorderingofparameterswhichisbrie ydiscussednextbeforeproceedingwithanexample.Parameternumbering.Whenusingtheconpatandfixedarguments,completeparam-etervectorsshouldbesupplied,i.e.,thesevectorsshouldhavelengthequaltothenumberofparametersofthemodel,whichcanbeobtainedbycallingnpar(object).Notethatthisisnotthesameasthedegreesoffreedomusede.g.,inthelogLikfunctionbecausenparalsocountsthebaselinecategoryzeroesfromthemultinomiallogisticmodels.Parametersarenumberedinthefollowingorder:1.thepriormodelparameters;2.theparametersforthetransitionmodels;3.theresponsemodelparametersperstate(andsubsequentlyperresponseinthecaseofmultivariatetimeseries).Toseetheorderingofparametersusethefollowing:R�setpars(mod,value=1:npar(mod))Toseewhichparametersarexed(bydefaultonlybaselineparametersinthemultinomiallogisticmodelsforthetransitionmodelsandtheinitialstateprobabilitiesmodel):R�setpars(mod,getpars(mod,which="fixed"))Whenttingconstraintsitisusefultohavegoodstartingvaluesfortheparametersandhencewersttthefollowingmodelwithoutconstraints: 12depmixS4:AnRPackageforHiddenMarkovModelsR�trst-c(0.9,0.1,0,0,0.1,0.9,0,0)Rကmod-depmix(list(rt~1,corr~1),data=speed,transition=~Pacc,+nstates=2,family=list(gaussian(),multinomial("identity")),+trstart=trst,instart=c(0.99,0.01))Rကfm1-fit(mod,verbose=FALSE,emc=em.control(rand=FALSE))convergedatiteration23withlogLik:-255.5Afterthis,weusethettedvaluesfromthismodeltoconstraintheregressioncoecientsonthetransitionmatrix(parametersnumber6and10):Rကpars-c(unlist(getpars(fm1)))Rကpars[6]-pars[10]-11Rကpars[1]-0Rကpars[2]-1Rကpars[13]-pars[14]-0.5Rကfm1-setpars(mod,pars)Rကconpat-c(0,0,rep(c(0,1),4),1,1,0,0,1,1,1,1)Rကconpat[6]-conpat[10]-2Rကfm2-fit(fm1,equal=conpat)Usingsummaryonthettedmodelshowsthattheregressioncoecientsarenowestimatedatthesamevalueof12.66.Settingelements13and14ofconpattozeroresultedinxingthoseparametersattheirstartingvaluesof0.5.Thismeansthatstate1cannowbeinterpretedasa'guessing'statewhichisassociatedwithcomparativelyfastresponses.Similarlyforelements1and2,resultinginxedinitialprobabilities.Thefunctionllratiocomputesthelikelihoodratio2-statisticandtheassociatedp-valuewithappropriatedegreesoffreedomfortestingthetenabilityofconstraints;Giudici,Ryden,andVandekerkhove(2000)showsthatthe`standardasymptotictheoryoflikelihood-ratiotestsisvalid'inhiddenMarkovmodels.(DannemannandHolzmann2007)discussesextensiontonon-standardsituations,e.g.fortestingparametersontheboundary.Notethatthesearguments(i.e.,conpatandconrows)providethepossibilityforarbitraryconstraints,alsobetween,e.g.,amultinomialregressioncoecientforthetransitionmatrixandthemeanofaGaussianresponsemodel.Whethersuchconstraintsmakesenseishencetheresponsibilityoftheuser.3.6.AddingcovariatesonthepriorprobabilitiesToillustratetheuseofcovariatesonthepriorprobabilitieswehaveincludedanotherdatasetwithdepmixS4.Thebalancedataconsistsof4binaryitems(correct-incorrect)onabalancescaletask(Siegler1981).ThedataformasubsetofthedatapublishedinJansenandvanderMaas(2002).Beforespecifyingspecifyingamodelforthesedata,webrie ydescribethem.ThebalancescaletaskisafamoustaskfortestingcognitivestrategiesdevelopedbyJeanPiaget(seeSiegler1981).Figure2providesanexampleofabalancescaleitem.Participants'taskistosaytowhichsidethebalancewilltipwhenreleased,oralternatively,whetheritwillstayinbalance.TheitemshowninFigure2isaso-calleddistanceitem:thenumberofweightsplacedoneachsideisequal,andonlythedistanceoftheweightstothefulcrumdiersbetweeneachside. IngmarVisser,MaartenSpeekenbrink13 Figure2:Balancescaleitem;thisisadistanceitem(seethetextfordetails).Childreninthelowergradesofprimaryschoolareknowntoignorethedistancedimension,andbasetheiransweronlyonthenumberofweightsoneachside.Hence,theywouldtypicallyprovidethewronganswertothesedistanceitems.Slightlyolderchildrendotakedistanceintoaccountwhenrespondingtobalancescaleitems,buttheyonlydosowhenthenumberofweightsisequaloneachside.ThesetwostrategiesthatchildrenemployareknownasRuleIandRuleII.Otherstrategiescanbeteasedapartbyadministeringdierentitems.Thebalancedatasetthatweanalysehereconsistsof4distanceitemsonabalancescaletaskadministeredto779participantsrangingfrom5to19yearsofage.Thefullsetofitemsconsistedof25items;otheritemsinthetestareusedtodetectotherstrategiesthatchildrenandyoungadultsemployinsolvingbalancescaleitems(seeJansenandvanderMaas2002,fordetails).Inthefollowingmodel,ageisincludedascovariateonclassmembershiptotestwhether,withage,childrenapplymorecomplexrulesinsolvingbalancescaleitems.Similarlytothetransitionmatrix,covariatesonthepriorprobabilitiesofthelatentstates(orclassesinthiscase),aredenedbyusingaone-sidedformulaprior=~age:R�data("balance")R�set.seed(1)R�mod-mix(list(d1~1,d2~1,d3~1,d4~1),data=balance,+nstates=3,family=list(multinomial("identity"),+multinomial("identity"),multinomial("identity"),+multinomial("identity")),respstart=runif(24),prior=~age,+initdata=balance)Rကfm-fit(mod,verbose=FALSE,emc=em.control(rand=FALSE))convergedatiteration77withlogLik:-917.5RကfmConvergenceinfo:Loglikelihoodconvergedtowithintol.(relativechange)'logLik.'-917.5(df=16)AIC:1867BIC:1942Noteherethatwedeneamixmodelinsteadofadepmixmodelasthesedataformindepen-dentobservations.Moreformally,depmixmodelsextendtheclassof`mix'modelsbyadding 14depmixS4:AnRPackageforHiddenMarkovModelstransitionmodels.Asforttingmixmodels:ascanbeseeninEquation9,theEMalgorithmcanbeappliedbysimplydroppingthesecondsummandcontainingthetransitionparameters,andthisisimplementedassuchintheEMalgorithmsindepmixS4.Ascanbeseenfromtheprintofthettedmodelabove,theBICforthismodelequals1941.6.Thesimilarlydened2-classmodelforthesedatahasaBICof1969.2,andthe4-classmodelhasBICequalto1950.4.Hence,the3-classseemstobeadequatefordescribingthesedata.Thesummaryofthefittedmodelgivesthefollowing(onlythepriormodelisshownhere):R�summary(fm,which="prior")MixtureprobabilitiesmodelModeloftypemultinomial(mlogit),formula:~ageCoefficients:St1St2St3(Intercept)06.39571.7548age0-0.6763-0.2905Probalitiesatzerovaluesofthecovariates.0.001650.98880.009541Theinterceptvaluesofthemultinomiallogisticparametersprovidethepriorprobabilitiesofeachclasswhenthecovariatehasvaluezero(notethatinthiscasethevaluezerodoesnothavemuchsignicance,scalingand/orcenteringofcovariatesmaybeusefulinsuchcases).Thesummaryfunctionprintsthesevalues.Ascanbeseenfromthosevalues,atagezero,thepriorprobabilityisoverwhelminglyatthesecondclass.Inspectionoftheresponseparametersrevealsthatclass2isassociatedwithincorrectresponding,whereasclass1isassociatedwithcorrectresponding;class3isanintermediateclasswithguessingbehavior.Figure3depictsthepriorclassprobabilitiesasfunctionofagebasedonthettedparameters.AscanbeseenfromFigure3,atyoungerageschildrenpredominantlyapplyRuleI,thewrongstrategyfortheseitems.Accordingtothemodel,approximately90%ofchildrenatage5applyRuleI.Theremaining10%areevenlysplitamongthe2otherclasses.Atage19,almostallparticipantsbelongtoclass1.Interestingly,priorprobabilityofthe'guess'class2,rstincreaseswithage,andthendecreasesagain.Thissuggeststhatchildrenpassthroughaphaseinwhichtheyareuncertainorpossiblyswitchbetweenapplyingdierentstrategies.4.ExtendingdepmixS4ThedepmixS4packagewasdesignedwiththeaimofmakingitrelativelyeasytoaddnewresponsedistributions(aswellaspossiblynewpriorandtransitionmodels).Tomakethispossible,theEMroutinesimplycallsthefitmethodsoftheseparateresponsemodelswithoutneedingaccesstotheinternalworkingsoftheseroutines.Referringtoequation9,theEMalgorithmcallsseparatetfunctionsforeachpartofthemodel,thepriorprobabilitymodel,thetransitionmodels,andtheresponsemodels.Asaconsequence,addinguser-speciedresponsemodelsisstraightforward.ThecurrentlyimplementeddistributionsarelistedinTable1.User-deneddistributionsshouldextendthe`response'classandhavethefollowingslots: IngmarVisser,MaartenSpeekenbrink15 Figure3:Classprobabilitiesasafunctionofage. packagefamilylink statsbinomiallogit,probit,cauchit,log,cloglogstatsgaussianidentity,log,inversestatsGammainverse,identity,logstatspoissonlog,identity,sqrtdepmixS4multinomiallogit,identity(nocovariatesallowed)depmixS4multivariatenormalidentity(onlyavailablethroughmakeDepmix)depmixS4ex-gaussidentity(onlyavailablethroughmakeDepmixasexample) Table1:ResponsedistributionavailableindepmixS4. 16depmixS4:AnRPackageforHiddenMarkovModels1.y:Theresponsevariable.2.x:Thedesignmatrix,possiblyonlyanintercept.3.parameters:Anamedlistwiththecoecientsandpossiblyotherparameters(e.g.,thestandarddeviationinthenormalresponsemodel).4.fixed:Avectoroflogicalsindicatingwhetherparametersarexed.5.npar:Numericalindicatingthenumberofparametersofthemodel.Inthespeeddataexample,itmaybemoreappropriatetomodeltheresponsetimeswithanexgausratherthanaGaussiandistribution.Todoso,werstdenean`exgaus'classextendingthe`response'class:R�setClass("exgaus",contains="response")Theso-denedclassnowneedsanumberofmethods:1.constructor:Functiontocreateinstancesoftheclasswithstartingvalues.2.show:Toprintthemodeltotheterminal.3.dens:Thefunctionthatcomputesthedensityoftheresponses.4.getparsandsetpars:Togetandsetparameters.5.predict:Togeneratepredictedvalues.6.fit:Functiontotthemodelusingposteriorweights(usedbytheEMalgorithm).Onlytheconstructorandthefitmethodsareprovidedhere;thecompletecodecanbefoundinthehelpleofthemakeDepmixfunction.Theexamplewiththeexgausdistributionusesthegamlssandgamlss.distpackages(RigbyandStasinopoulos2005;StasinopoulosandRigby2007;Stasinopoulos,Rigby,andAkantziliotou2009;Stasinopoulos,Rigby,Akantzil-iotou,Heller,Ospina,andMotpan2010)forcomputingthedensityandforfittingtheparameters.Theconstructormethodreturnanobjectofclass`exgaus',andisdenedasfollows:R�library("gamlss")R�library("gamlss.dist")R�setGeneric("exgaus",function(y,pstart=NULL,fixed=NULL,...)+standardGeneric("exgaus"))R�setMethod("exgaus",+signature(y="ANY"),+function(y,pstart=NULL,fixed=NULL,...){+y-matrix(y,length(y))+x-matrix(1)+parameters-list()+npar-3+if(is.null(fixed))fixed-as.logical(rep(0,npar)) IngmarVisser,MaartenSpeekenbrink17+if(!is.null(pstart)){+if(length(pstart)!=npar)stop("lengthof'pstart'mustbe",npar)+parameters$mu-pstart[1]+parameters$sigma-log(pstart[2])+parameters$nu-log(pstart[3])+}+mod-new("exgaus",parameters=parameters,fixed=fixed,+x=x,y=y,npar=npar)+mod+}+)Thefitmethodisdenedasfollows:RကsetMethod("fit","exgaus",+function(object,w){+if(missing(w))w-NULL+y-object@y+fit-gamlss(y~1,weights=w,family=exGAUS(),+control=gamlss.control(n.cyc=100,trace=FALSE),+mu.start=object@parameters$mu,+sigma.start=exp(object@parameters$sigma),+nu.start=exp(object@parameters$nu))+pars-c(fit$mu.coefficients,fit$sigma.coefficients,+fit$nu.coefficients)+object-setpars(object,pars)+object+}+)Thefitmethoddenesagamlssmodelwithonlyanintercepttobeestimatedandthensetsthettedparametersbackintotheirrespectiveslotsinthe`exgaus'object.Seethehelpforgamlss.distrforinterpretationoftheseparameters.Afterdeningallthenecessarymethodsforthenewresponsemodel,wecannowdenethedependentmixturemodelusingthisresponsemodel.ThefunctionmakeDepmixisincludedindepmixS4tohavefullcontrolovermodelspecication,andweneedithere.Werstcreatealltheresponsemodelsthatweneedasadoublelist:RကrModels-list()RကrModels[[1]]-list()RကrModels[[1]][[1]]-exgaus(speed$rt,pstart=c(5,0.1,0.1))RကrModels[[1]][[2]]-GLMresponse(formula=corr~1,data=speed,+family=multinomial(),pstart=c(0.5,0.5))RကrModels[[2]]-list()RကrModels[[2]][[1]]-exgaus(speed$rt,pstart=c(6,0.1,0.1))RကrModels[[2]][[2]]-GLMresponse(formula=corr~1,data=speed,+family=multinomial(),pstart=c(0.1,0.9)) 18depmixS4:AnRPackageforHiddenMarkovModelsNext,wedenethetransitionandpriorprobabilitymodelsusingthetransInitfunction(whichproducesatransInitmodel,whichalsoextendsthe`response'class):R�trstart-c(0.9,0.1,0.1,0.9)Rကtransition-list()Rကtransition[[1]]-transInit(~Pacc,nst=2,data=speed,+pstart=c(0.9,0.1,0,0))Rကtransition[[2]]-transInit(~Pacc,nst=2,data=speed,+pstart=c(0.1,0.9,0,0))RကinMod-transInit(~1,ns=2,pstart=c(0.1,0.9),+data=data.frame(1))Finally,weputeverythingtogetherusingmakeDepmixandtthemodel:Rကmod-makeDepmix(response=rModels,transition=transition,+prior=inMod,homogeneous=FALSE)Rကfm-fit(mod,verbose=FALSE,emc=em.control(rand=FALSE))convergedatiteration43withlogLik:-232.3Usingsummarywillprintthettedparameters.NotethattheuseofmakeDepmixallowsthepossibilityof,say,ttingagaussianinonestateandanexgausdistributioninanotherstate.NotealsothataccordingtotheAICandBIC,themodelwiththeexgausdescribesthedatamuchbetterthanthesamemodelinwhichtheresponsetimesaremodelledasgaussian.5.ConclusionsandfutureworkdepmixS4providesa exibleframeworkforttingdependentmixturemodelsforalargevarietyofresponsedistributions.Itcanalsotlatentclassregressionandnitemixturemodels,althoughitshouldbenotedthatmorespecializedpackagesareavailableforthissuchas exmix(Leisch2004;GrunandLeisch2008).Thepackageisintendedformodelling(individual)timeseriesdatawiththeaimofcharacterizingthetransitionprocessesunderlyingthedata.Thepossibilitytousecovariatesonthetransitionmatrixgreatlyenhancesthe exibilityofthemodel.TheEMalgorithmusesaverygeneralinterfacethatallowseasyadditionofnewresponsemodels.Wearecurrentlyworkingonimplementingthegradientsforresponseandtransitionmod-elswithtwogoalsinmind.First,tospeedup(constrained)parameteroptimizationusingRdonlp2orRsolnp.Second,analyticgradientsareusefulincomputingtheHessianoftheestimatedparameterssoastoarriveatstandarderrorsforthose.Wearealsoplanningtoimplementgoodness-of-tstatistics(TitmanandSharples2008).AcknowledgmentsIngmarVisserwassupportedbyanECFramework6grant,project516542(NEST).MaartenSpeekenbrinkwassupportedbyESRCgrantRES-062-23-1511andtheESRCCentrefor IngmarVisser,MaartenSpeekenbrink19EconomicLearningandSocialEvolution(ELSE).HanvanderMaasprovidedthespeed-accuracydata(Dutilhetal.2011)andtherebynecessitatedimplementingmodelswithtime-dependentcovariates.BrendaJansenprovidedthebalancescaledataset(JansenandvanderMaas2002)whichwastheperfectopportunitytotestthecovariatesonthepriormodelparameters.Theexamplesinthehelplesusebothofthesedatasets.ReferencesAgrestiA(2002).CategoricalDataAnalysis.2ndedition.JohnWiley&Sons,Hoboken.BaumLE,PetrieT(1966).\StatisticalInferenceforProbabilisticFunctionsofFiniteStateMarkovChains."AnnalsofMathematicalStatistics,67,1554{40.CappeO,MoulinesE,RydenT(2005).InferenceinHiddenMarkovModels.Springer-Verlag,NewYork.ChungH,WallsT,ParkY(2007).\ALatentTransitionModelWithLogisticRegression."Psychometrika,72(3),413{435.DannemannJRN,HolzmannH(2007).\LikelihoodRatioTestingforHiddenMarkovModelsunderNon-StandardConditions."ScandinavianJournalofStatistics,25(2),309{321.DutilhG,WagenmakersEJ,VisserI,vanderMaasHLJ(2011).\APhaseTransitionModelfortheSpeed{AccuracyTrade{OinResponseTimeExperiments."CognitiveScience,35,211{250.Fruhwirth-SchnatterS(2006).FiniteMixtureandMarkovSwitchingModels.Springer-Verlag,NewYork.GhalanosA,TheulS(2010).Rsolnp:GeneralNon-LinearOptimizationUsingAugmentedLagrangeMultiplierMethod.Rpackageversion1.0-4,URLhttp://CRAN.R-project.org/package=Rsolnp.GhyselsE(1994).\OnthePeriodicStructureoftheBusinessCycle."JournalofBusinessandEconomicStatistics,12(3),289{298.GiudiciP,RydenT,VandekerkhoveP(2000).\Likelihood-RatioTestsforHiddenMarkovModels."Biometrics,56(3),742{747.GrunB,LeischF(2008).\FlexMixVersion2:FiniteMixtureswithConcomitantVariablesandVaryingandConstantParameters."JournalofStatisticalSoftware,28(4),1{35.URLhttp://www.jstatsoft.org/v28/i04/.JansenBRJ,vanderMaasHLJ(2002).\TheDevelopmentofChildren'sRuleUseontheBalanceScaleTask."JournalofExperimentalChildPsychology,81(4),383{416.KimCJ(1994).\DynamicLinearModelswithMarkov-Switching."JournalofEconometrics,60,1{22. 20depmixS4:AnRPackageforHiddenMarkovModelsKroghA(1998).\AnIntroductiontoHiddenMarkovModelsforBiologicalSequences."InSLSalzberg,DBSearls,SKasif(eds.),ComputationalMethodsinMolecularBiology,chapter4,pp.45{63.Elsevier,Amsterdam.LeischF(2004).\FlexMix:AGeneralFrameworkforFiniteMixtureModelsandLatentClassRegressioninR."JournalofStatisticalSoftware,11(8),1{18.URLhttp://www.jstatsoft.org/v11/i08/.LerouxBG,PutermanML(1992).\Maximum-Penalized-LikelihoodEstimationforIndepen-dentandMarkov-DependentMixtureModels."Biometrics,48,545{548.LystigTC,HughesJP(2002).\ExactComputationoftheObservedInformationMatrixforHiddenMarkovModels."JournalofComputationalandGraphicalStatistics,11(3),678{689.RabinerLR(1989).\ATutorialonHiddenMarkovModelsandSelectedApplicationsinSpeechRecognition."ProceedingsofIEEE,77(2),267{295.RainerG,MillerEK(2000).\NeuralEnsembleStatesinPrefrontalCortexIdentiedUsingaHiddenMarkovModelwithaModiedEMAlgorithm."Neurocomputing,32{33,961{966.RatcliR(1978).\ATheoryofMemoryRetrieval."PsychologicalReview,85,59{108.RDevelopmentCoreTeam(2010).R:ALanguageandEnvironmentforStatisticalComputing.RFoundationforStatisticalComputing,Vienna,Austria.ISBN3-900051-07-0,URLhttp://www.R-project.org/.RigbyRA,StasinopoulosDM(2005).\GeneralizedAdditiveModelsforLocation,ScaleandShape."AppliedStatistics,54,507{554.SchmittmannVD,VisserI,RaijmakersMEJ(2006).\MultipleLearningModesintheDe-velopmentofRule-basedCategory-learningTaskPerformance."Neuropsychologia,44(11),2079{2091.SieglerRS(1981).DevelopmentalSequenceswithinandbetweenConcepts.Number46inMonographsoftheSocietyforResearchinChildDevelopment.SRCD.SpellucciP(2002).donlp2UsersGuide.TUDarmstadt.URLhttp://www.mathematik.tu-darmstadt.de/fbereiche/numerik/staff/spellucci/DONLP2/.StasinopoulosDM,RigbyBA(2007).\GeneralizedAdditiveModelsforLocationScaleandShape(GAMLSS)inR."JournalofStatisticalSoftware,23(7),1{46.URLhttp://www.jstatsoft.org/v23/i07.StasinopoulosDM,RigbyBA,AkantziliotouC(2009).gamlss:GeneralizedAdditiveModelsforLocationScaleandShape.Rpackageversion2.0-0,URLhttp://CRAN.R-project.org/package=gamlss.StasinopoulosDM,RigbyBA,AkantziliotouC,HellerG,OspinaR,MotpanN(2010).gamlss.dist:DistributionstoBeUsedforGAMLSSModelling.Rpackageversion4.0-0,URLhttp://CRAN.R-project.org/package=gamlss.dist. IngmarVisser,MaartenSpeekenbrink21TamuraR(2009).Rdonlp2:AnRExtensionLibrarytoUsePeterSpelluci'sdonlp2fromR.Rpackageversion0.4,URLhttp://arumat.net/Rdonlp2/.TitmanAC,SharplesLD(2008).\AGeneralGoodness-of-tTestforMarkovandHiddenMarkovModels."StatisticsinMedicine,27,2177{2195.vandePolF,LangeheineR,JongWD(1996).Panmark3.PanelAnalysisUsingMarkovChains.ALatentClassAnalysisProgram.Voorburg.VanderMaasHLJ,MolenaarPCM(1992).\StagewiseCognitiveDevelopment:AnApplica-tionofCatastropheTheory."PsychologicalReview,99,395{417.VenablesWN,RipleyBD(2002).ModernAppliedStatisticswithS.4thedition.Springer-Verlag,NewYork.VermuntJK,MagidsonJ(2003).LatentGold3.0.Belmont,MA.URLhttp://www.statisticalinnovations.com/.VisserI,RaijmakersMEJ,VanderMaasHLJ(2009).\HiddenMarkovModelsforIndividualTimeSeries."InJValsiner,PCMMolenaar,MCDPLyra,NChaudhary(eds.),DynamicProcessMethodologyintheSocialandDevelopmentalSciences,chapter13,pp.269{289.Springer-Verlag,NewYork.VisserI,SpeekenbrinkM(2010).\depmixS4:AnRPackageforHiddenMarkovModels."JournalofStatisticalSoftware,36(7),1{21.URLhttp://www.jstatsoft.org/v36/i07/.WickensTD(1982).ModelsforBehavior:StochasticProcessesinPsychology.W.H.FreemanandCompany,SanFrancisco.YeY(1987).InteriorAlgorithmsforLinear,Quadratic,andLinearlyConstrainedNon-LinearProgramming.Ph.D.thesis,DepartmentofESS,StanfordUniversity.ZucchiniW,MacDonaldI(2009).HiddenMarkovModelsforTimeSeries:AnIntroductionUsingR.MonographsonStatisticsandAppliedProbability.CRCPress,BocaRaton.Aliation:IngmarVisserDepartmentofPsychologyUniversityofAmsterdamRoetersstraat151018WB,Amsterdam,TheNetherlandsE-mail:i.visser@uva.nlURL:http://www.ingmar.org/