### Code for conducting SEM in R - Analysis of Karels et al. data

### A big thank-you to Steve Dobson for providing the data!

### Uses package 'sem'

library(sem)

### import data
datum=read.csv(file.choose())
head(datum)

###Log transform the data (as Karels et al. did)
### Be sure to add '1' to all count data - can't log transform 0!
datum$NativeBirds=datum$NativeBirds+1
datum$TotalExtinctions=datum$TotalExtinctions+1
datum$Predators=datum$Predators+1
datumLN=log(datum)
head(datumLN)


### first step is to specify the path model, which is done with 'specify.model'
###Build the path diagram
###Arguments are path, variable name, starting value (NA lets R choose starting value)
modelPop=specifyModel()
	Area<>Area, AreaVar, NA ### variance in Area
	Isolation<>Isolation, IsoVar, NA ###variance in Isolation
	NativeBirds<>NativeBirds, BirdVar, NA ### error in No. of avifauna
	Predators<>Predators, PredVar, NA ### error in No. of predators
	TotalExtinctions<>TotalExtinctions, ExtVar, NA ### error in No. of Extinctions
	Area<>Isolation, AreaIso, NA ### correlation between area and isolation
	Area>Predators, AreaPred, NA ### link from Area to No. of Predators
	Area>NativeBirds, AreaBirds, NA ### link from Area to No. of avifauna
	Area>TotalExtinctions, AreaExt, NA ### link from Area to No. of Extinctions
	Isolation>NativeBirds, IsoBirds, NA ### link from Isolation to No. of avifauna
	Isolation>TotalExtinctions, IsoExt, NA ### link from isloation to No. of Extinctions
	Isolation>Predators, IsoPred, NA ### link from isolation to No. of Predators
	NativeBirds>TotalExtinctions, BirdsExt, NA ### link from No. of avifauna to No. of Extinctions
	Predators>TotalExtinctions, PredExt, NA ### link from No. of Predators to No. of Extinctions



### Run the SEM

results=sem(modelPop,data=datumLN)


### summary provides goodness-of-fit tests
### parameter estimates for variances/errors are standard deviation estimates of residuals
### parameter estimates for paths are betas of coefficients in regression equations
### compare estiamtes to truth (in excel data-sheet)
### NOTE - THESE ARE NOT STANDARDIZED PARAMETER COEFFICIENTS 
summary(results)


### function std.coef provides standardized path coefficiets
help(stdCoef)
stdCoef(results)

#                    Std. Estimate                                       
#AreaVar     AreaVar    1.00000000                         Area <--> Area
#IsoVar       IsoVar    1.00000000               Isolation <--> Isolation
#BirdVar     BirdVar    0.41289643           NativeBirds <--> NativeBirds
#PredVar     PredVar    0.58382583               Predators <--> Predators
#ExtVar       ExtVar    0.70615842 TotalExtinctions <--> TotalExtinctions
#AreaIso     AreaIso   -0.27777267                    Isolation <--> Area
#AreaPred   AreaPred    0.67135044                    Predators <--- Area
#AreaBirds AreaBirds    0.63469743                  NativeBirds <--- Area
#AreaExt     AreaExt    0.08455985             TotalExtinctions <--- Area
#IsoBirds   IsoBirds   -0.28775120             NativeBirds <--- Isolation
#IsoExt       IsoExt    0.28103045        TotalExtinctions <--- Isolation
#IsoPred     IsoPred    0.17103616               Predators <--- Isolation
#BirdsExt   BirdsExt    0.36630441      TotalExtinctions <--- NativeBirds
#PredExt     PredExt    0.22452040        TotalExtinctions <--- Predators
### Path coefficients for variances/errors are % of data not explained by path model (1-r^2)
### Path coefficients for paths are correlation coefficients (r) AFTER accounting for other effects	
### In essence, path coefficients for paths are relationship between variables when all other variables are fixed

### results should match numbers in Karels et al. pretty close (within about 0.02)


### What does this mean?

#ExtVar = 0.706; so paths explain about 30% of variation in total extinctions (70% explained by error)
#paths explain about 59% of variation in number of birds (1-BirdVar)
#paths explain about 42% of variation in number of predators (1-PredVar) 

### after statistical controlling for effects of isolation, area, and number of birds,
### predators had a moderate effect on extinctions (correlation - r; PredExt=0.22)