## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(pwSEM)

## -----------------------------------------------------------------------------
temp<-data.frame(X=rnorm(6),Y=rnorm(6),Z=rnorm(6),YEAR=c(2003,2003,2004,2004,2005,2006),
                 POP=c("A","B","A","C","B","C"))
head(temp)

## ----eval=FALSE---------------------------------------------------------------
# my.list<-list(mgcv::gam(X1~1,data=sim_normal.no.nesting,family=gaussian),
#          mgcv::gam(X2~X1,data=sim_normal.no.nesting,family=gaussian),
#          mgcv::gam(X3~X2,data=sim_normal.no.nesting,family=gaussian),
#          mgcv::gam(X4~X3,data=sim_normal.no.nesting,family=gaussian))

## ----eval=FALSE---------------------------------------------------------------
# out<-pwSEM(sem.functions=my.list,
#     marginalized.latents=list(X4~~X2), data=sim_normal.no.nesting,use.permutations = TRUE,n.perms=5000,do.smooth=FALSE,all.grouping.vars = NULL)

## ----eval=FALSE---------------------------------------------------------------
# summary(out,structural.equations=TRUE)

## ----eval=FALSE---------------------------------------------------------------
# my.list<-list(mgcv::gam(X1~1,data=sim_poisson.no.nesting,family=gaussian),
#          mgcv::gam(X2~X1,data=sim_poisson.no.nesting,family=poisson),
#          mgcv::gam(X3~X2,data=sim_poisson.no.nesting,family=poisson),
#          mgcv::gam(X4~X3,data=sim_poisson.no.nesting,family=poisson))
# out<-pwSEM(sem.functions=my.list,
#            marginalized.latents=list(X4~~X2),   data=sim_poisson.no.nesting,use.permutations = TRUE,n.perms=10000)
# summary(out,structural.equations=TRUE)

## -----------------------------------------------------------------------------
my.list<-list(gamm4::gamm4(X1~1,random=~(1|group),data=sim_normal.with.nesting,family=gaussian),
         gamm4::gamm4(X2~X1,random=~(1|group),data=sim_normal.with.nesting,family=gaussian),
         gamm4::gamm4(X3~X2,random=~(1|group),data=sim_normal.with.nesting,family=gaussian),
         gamm4::gamm4(X4~X3,random=~(1|group),data=sim_normal.with.nesting,family=gaussian))

out<-pwSEM(sem.functions=my.list,
           marginalized.latents=list(X4~~X2),
   data=sim_normal.with.nesting,use.permutations = TRUE,
   do.smooth=TRUE,all.grouping.vars=c("group"))
summary(out,structural.equations=TRUE)

## -----------------------------------------------------------------------------
library(ggm)
full.dag<-DAG(X2~X1+L1,X3~X2,X4~X3+L1,L2~X4+X5)
drawGraph(full.dag)
equivalent.dag<-DAG.to.MAG.in.pwSEM(full.DAG=full.dag,
                latents=c("L1","L2"),conditioning.latents=
                  c("L2"))
drawGraph(equivalent.dag)


## -----------------------------------------------------------------------------
library(ggm)
mag<-makeMG(dg=DAG(X2~X1,X3~X2,X4~X3),bg=UG(~X2*X4))
drawGraph(mag)
equivalent.dag<-MAG.to.DAG.in.pwSEM(MAG=mag,
       marginalized.latents=list(X2~~X4),conditioned.latents= NULL)
drawGraph(equivalent.dag)

## -----------------------------------------------------------------------------
R1<-residuals(mgcv::gam(X1~X2,data=sim_poisson.no.nesting,family=gaussian))
R2<-residuals(mgcv::gam(X3~X2,data=sim_poisson.no.nesting,family=poisson))
generalized.covariance(R1,R2)

## -----------------------------------------------------------------------------
R1<-residuals(mgcv::gam(X1~X2,data=sim_poisson.no.nesting,family=gaussian))
R2<-residuals(mgcv::gam(X3~X2,data=sim_poisson.no.nesting,family=poisson))
perm.generalized.covariance(R1,R2,nperm=5000)

## -----------------------------------------------------------------------------
# DAG: X1->X2->X3->X4 and X2<->X4
my.list<-list(mgcv::gam(X1~1,data=sim_poisson.no.nesting,family=gaussian),
         mgcv::gam(X2~X1,data=sim_poisson.no.nesting,family=poisson),
          mgcv::gam(X3~X2,data=sim_poisson.no.nesting,family=poisson),
          mgcv::gam(X4~X3,data=sim_poisson.no.nesting,family=poisson))
 out<-pwSEM(sem.functions=my.list,marginalized.latents=list(X4~~X2),
           data=sim_poisson.no.nesting,use.permutations = FALSE)
view.paths(from="X1",to="X4",sem.functions=out$sem.functions,data=
sim_poisson.no.nesting,scale="response",dag=out$causal.graph)

## -----------------------------------------------------------------------------
library(ggm)
mag<-makeMG(dg=DAG(X2~X1,X3~X2,X4~X3),bg=UG(~X2*X4))
mag
drawGraph(mag)
basiSet.MAG(cgraph=mag)

## -----------------------------------------------------------------------------
mag<-makeMG(dg=DAG(X3~X1+X2,X4~X3),bg=UG(~X1*X2))
drawGraph(mag)
set.seed(11)
L<-rnorm(200)
X1<-0.5*L+rnorm(200,0,sqrt(1-0.5^2))
X2<-0.5*L+rnorm(200,0,sqrt(1-0.5^2))
X3<-0.5*X1+0.5*X2+rnorm(200,0,sqrt(1-2*0.5^2))
X4<-0.5*X3+rnorm(200,0,sqrt(1-0.5^2))
my.dat<-data.frame(X1,X2,X3,X4)
my.list<-list(mgcv::gam(X1~1,data=my.dat,family=gaussian),
        mgcv::gam(X2~1,data=my.dat,family=gaussian),
        mgcv::gam(X3~X1+X2,data=my.dat,family=gaussian),
        mgcv::gam(X4~X3,data=my.dat,family=gaussian))
get.AIC(sem.model=my.list,MAG=mag,data=my.dat)

## -----------------------------------------------------------------------------
MCX2(model.df=4,n.obs=15,model.chi.square=9.89)

## ----eval=FALSE---------------------------------------------------------------
# nesting.structure<-list(XR=c("year","nest"),XM=c("year","nest"),
#              XH=c("year","nest"),XP=c("year","nest"),
#              XF=c("year","nest"))
# CI.algorithm(dat=nested_data[,-3],family=data.frame(XR="binomial",
#             XP="poisson"),nesting=nesting.structure,alpha.reject=0.3,
#             smooth=T)

## -----------------------------------------------------------------------------
nesting.structure<-list(XR=c("year","nest"),XM=c("year","nest"),
             XH=c("year","nest"),XP=c("year","nest"),
             XF=c("year","nest"))
con.edges<-"
XF|XM
XH->XR
"

CI.algorithm(dat=nested_data[,-3],family=data.frame(XR="binomial",
            XP="poisson"),nesting=nesting.structure,alpha.reject=0.3,
            constrained.edges=con.edges,smooth=T)

## ----eval=FALSE---------------------------------------------------------------
# CI.algorithm(sim_tetrads)

## -----------------------------------------------------------------------------
vanishing.tetrads(dat=sim_tetrads,sig=0.05)

## -----------------------------------------------------------------------------
vanishing.tetrads(dat=sim_tetrads,sig=0.05,bootstrap=TRUE,B=1000)