* Experiment 1

version 16
clear all
macro drop _all
set linesize 80
set more off
set scheme s1manual
set graphics off /* create graphs without them popping up on screen */

//----------- open data file
use data/Expt1_2023_forAnalysis_withSyll.dta, clear
di _N

// how many trials ?
preserve
bysort subject trial: keep if _n == 1
di _N
restore


//# ----- INITIATION
//# -- Descriptive stats
tabstat keytime if keytime > 50 & keytime < 1500 & keyNum == 1, by(wordtype)  statistics(n mean sem sd min max) columns(statistics)
tabstat keytime if keytime > 50 & keytime < 1500 & keyNum == 2, by(wordtype)  statistics(n mean sem sd min max) columns(statistics)


*MODEL1: cpm crump. k1 vs k2  set wordtype wlen wfreq
* model 1
mixed keytime i.wordtypeN##i.keyNum i.keyNum##c.neWord_SLlg10wf i.keyNum##c.wordlen || _all:R.subject || item: if  isControl != 1 & corLD == 1 & keytime > 50 & keytime < 1500 & inlist(keyNum, 1,2) , mle nogroup nolrtest nofvlabel noretable
contrast  i.wordtypeN##i.keyNum
contrast i.keyNum#c.neWord_SLlg10wf
margins , dydx(neWord_SLlg10wf)
contrast i.keyNum#c.wordlen
margins i.keyNum, dydx(c.wordlen)



*MODEL 2 and 3:  ROLE OF CONSTITUENTS INFO ON K1 C AND P ONLY
*Remember that c1len c2len and wordlen can't all be in same model due to collinearity

** SEPARATE MODELS FOR C1 AND C2 VARIABLES
*with just c1 variables
mixed keytime ///
	i.wordtypeN##c.neC1SLlg10wf i.wordtypeN##c.c1len ///
	|| _all:R.subject || item: if keytime > 50 & keytime < 1500 & inlist(keyNum, 1) & isControl != 1 , mle

contrast i.wordtypeN#c.neC1SLlg10wf
contrast  i.wordtypeN#c.c1len
margins i.wordtypeN, dydx(c.c1len)


*figure  "Predicted initiation latencies in Experiment 1 as a function of wordtype and length of the first constituent"
margin wordtypeN, at(c1len = (3(1)6))
marginsplot,  ytitle ("Predicted latencies (ms)") xtitle("Length of first constituent") title("") legend(subtitle("Wordtype") rows(2) ring(0) pos(11)) name(E1_i_c1len, replace) aspectratio(.75)
graph save  "output/figures/fig_E1_init_c1len" , replace
graph export "output/figures/fig_E1_init_c1len.png", replace*/

* with c2 variables
mixed keytime ///
	i.wordtypeN##c.nec2SLlg10wf i.wordtypeN##c.c2len ///
	|| _all:R.subject || item: if keytime > 50 & keytime < 1500 & inlist(keyNum, 1) & isControl != 1 , mle nogroup nolrtest nofvlabel noretable

contrast i.wordtypeN#c.nec2SLlg10wf
margins , dydx(c.nec2SLlg10wf)
contrast  i.wordtypeN#c.c2len
margins , dydx(c.c2len)




//# ----- BOUNDARY
//# -- Descriptive stats
tabstat keytime if keytime > 50 & keytime < 1000 & dist == -1, by(wordtype)  statistics(n mean sem sd min max) columns(statistics)
tabstat keytime if keytime > 50 & keytime < 1000 &  dist == 0, by(wordtype)  statistics(n mean sem sd min max) columns(statistics)

mixed keytime  i.distNum##i.wordtypeN || _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(dist,-1,0)  , mle nogroup nolrtest nofvlabel noretable
contrast i.distNum##i.wordtypeN
contrast i.distNum@i.wordtypeN

*Is the Boundary Effect Influenced by Word and Orthographic Properties?
* INITIAL MODEL c p include wordlevel variables and bigram (mentioned but not reported in detail)
mixed keytime  i.distNum##i.wordtypeN##c.neWord_SLlg10w ///
	i.distNum##i.wordtypeN##c.bgJonesMewhortboundaryNE ///
	|| _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(dist,-1,0) & wordtypeN != 2 , mle nogroup nolrtest nofvlabel noretable
contrast i.distNum##i.wordtypeN#c.neWord_SLlg10w
margins, dydx(neWord_SLlg10w) /* marginal */
* no evidence that wordfreq interacts so remove interaction

* FINAL MODEL (reported)
mixed keytime c.neWord_SLlg10w ///
	i.distNum##i.wordtypeN##c.bgJonesMewhortboundaryNE ///
	|| _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(dist,-1,0) & wordtypeN != 2 , mle nogroup nolrtest nofvlabel noretable

margins, dydx(neWord_SLlg10w) /*sig*/

* -1 of p not affected by bigram
contrast i.distNum#i.wordtypeN#c.bgJonesMewhortboundaryNE
contrast i.distNum#c.bgJonesMewhortboundaryNE@wordtypeN
margins i.distNum#3.wordtypeN, dydx(c.bgJonesMewhortboundaryNE)
* for compounds
margins 1.wordtypeN, dydx(c.bgJonesMewhortboundaryNE)

*Do properties of the (pseudo) constituent influence processing at the boundary region, and if so does morphological structure influence this?

*C1 Properties
* --- MODEL -1,0 cp constituent properties of c1
mixed keytime c.neWord_SLlg10w i.distNum##i.wordtypeN##c.bgJonesMewhortboundaryNE ///
	i.distNum##i.wordtypeN##c.neC1SLlg10wf i.distNum##i.wordtypeN##c.c1len || _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(dist,-1,0) , mle nogroup nolrtest nofvlabel noretable


* length
contrast i.distNum#i.wordtypeN#c.c1len

* do 0 and -1 slopes differ for p and c.
* see p. 48 https://www.stata.com/manuals13/rcontrast.pdf

contrast r.distNum#c.c1len@wordtypeN
* do slope differs from 0 for pos 0 and -1 for c
margins distNum, dydx(c1len) subpop(if wordtypeN==1)

* do slope differs from 0 for pos 0 and -1 for p
margins distNum, dydx(c1len) subpop(if wordtypeN==3)

* no effect of c1freq
contrast i.distNum#i.wordtypeN#c.neC1SLlg10wf
contrast i.wordtypeN#c.neC1SLlg10wf
contrast i.distNum#c.neC1SLlg10wf
margins, dydx(neC1SLlg10wf)

*figure
quietly margins distNum#wordtypeN, at(c1len=(3/6))
marginsplot,  ytitle ("Predicted latencies (ms)") xtitle("Length of first constituent") title("") by(wordtypeN) legend( order(3 "end of C1" 4 "start of C2") rows(2) ring(0) pos(5)) byopts(title(" ")) subtitle(, bcolor(none) ) ylabel(150(50)350)  name(E1_b_c1len, replace) aspectratio(.75)

graph save  "output/figures/fig_E1_b_c1len", replace
graph export  "output/figures/fig_E1_b_c1len.png", replace

* no effect of frequency so not useful to plot
quietly margins distNum#wordtypeN, at(neC1SLlg10wf=(0.3/5))
marginsplot,  by( wordtypeN ) xtitle(Frequency of the first constituent) ytitle(Predicted latencies (ms)) byopts(title(" ") legend(off)) subtitle(, bcolor(none) ) legend(order(3 "end of C1" 4 "start of C2" )) ylabel(150(50)350) name(E1_b_c1freq, replace) aspectratio(.75)

* combined and save as stata & png formats
grc1leg2 E1_b_c1len E1_b_c1freq, cols(1) ///
	imargin(0 0 0 0) ycommon legendfrom(E1_b_c1len)  xsize(5) ysize(5)  pos(5)
graph save  "output/figures/fig_E1_b_c1", replace
graph export  "output/figures/fig_E1_b_c1.png", replace


*C2 Properties

* --- MODEL -1,0 cp constituent properties of c2

mixed keytime c.neWord_SLlg10w i.distNum##i.wordtypeN##c.bgJonesMewhortboundaryNE ///
	i.distNum##i.wordtypeN##c.nec2SLlg10wf i.distNum##i.wordtypeN##c.c2len || _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(dist,-1,0) , mle nogroup nolrtest nofvlabel noretable

* c2len
contrast i.distNum#i.wordtypeN#c.c2len
contrast i.distNum#c.c2len@wordtypeN

* do slopes differ from 0 for p
margins distNum, dydx(c2len) subpop(if wordtypeN==3)
* do slopes differ from 0 for c
margins , dydx(c2len) subpop(if wordtypeN==1)
* do 0 and -1 slopes differ for p and c
contrast r.distNum#c.c2len@wordtypeN

*c2freq
contrast i.distNum#i.wordtypeN#c.nec2SLlg10wf
contrast i.wordtypeN#c.nec2SLlg10wf
contrast i.distNum#c.nec2SLlg10wf
margins distNum, dydx(nec2SLlg10wf)
margins , dydx(nec2SLlg10wf)


quietly margins distNum#wordtypeN, at(c2len=(3/8))
marginsplot,  ytitle ("Predicted latencies (ms)") xtitle("Length of second constituent") title("") by(wordtypeN) legend( order(3 "end of C1" 4 "start of C2") rows(2) ring(0) pos(5)) byopts(title(" ")) subtitle(, bcolor(none) ) ylabel(150(50)350)  name(E1_b_c2len, replace)

quietly margins distNum, at(nec2SLlg10wf=(0.3/5))
marginsplot,   xtitle(Frequency of the second constituent) ytitle(Predicted latencies (ms)) title(" ") byopts(title(" ") legend(off)) subtitle(, bcolor(none) ) legend(order(3 "end of C1" 4 "start of C2" )) ylabel(150(50)350) name(E1_b_c2freq, replace)  aspectratio(.75)

grc1leg2 E1_b_c2len E1_b_c2freq, rows(2)  ///
	imargin(0 0 0 0) ycommon legendfrom(E1_b_c2len) ysize(5)   pos(5)
graph save  "output/figures/Fig2fig_E1_b_c2", replace
graph export  "output/figures/Fig2fig_E1_b_c2.png", replace


//# Is the morpheme boundary effect merely a syllable effect?
cd "/Users/christinagagne/Dropbox/Analyses/2023_JEPLMC__analysis"
use data/Expt1_2023_forAnalysis_withSyll.dta, clear

*model 1 We focused on the compounds and pseudo-compounds where the morpheme boundary coincided with a syllable and fit a model with position (last letter of the first constituent vs. first letter of the second constituent) and wordtype (compound and pseudo-compound) as fixed factors and subjects and items as crossed-random effects.

mixed keytime i.distNum##i.wordtypeN || _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(dist,-1,0) & wordtypeN != 2 & atSyll == 1, mle nogroup nolrtest nofvlabel noretable
contrast i.distNum##i.wordtypeN
contrast i.distNum@i.wordtypeN

* same as interaction in this case because 2x2
contrast {distNum -1 1}#{wordtypeN -1 1}

quietly margins distNum#wordtypeN
marginsplot, noci ytitle ("Predicted latencies (ms)") xtitle("Position relative to the first letter of the second constituent") title("") legend( order(1 "compound" 2 "pseudo-compound") rows(2) ring(0) pos(11)) byopts(title(" ")) subtitle(, bcolor(none) )
graph save  "output/figures/fig_E1_b_atsyllone", replace
graph export  "output/figures/fig_E1_b_atsyllone.png", replace



*model 2 To further investigate the potential influence of morphemes vs. syllables, we examined whether there is any additional influence of a syllable boundary corresponding to a pseudo-morpheme boundary or not.

preserve
keep if wordtypeN == 3
drop if firstsyll == 0 & firstsyll_pre == 0
gen distsyll = -1 if dist == -1 & firstsyll == 1
replace distsyll = 0 if dist == 0 & firstsyll == 1
*-M+S items shift up one position
replace distsyll = -1 if dist == -2 & firstsyll == 0
replace distsyll = 0 if dist == -1 & firstsyll == 0
egen distsyllNum = group(distsyll), label(distsyllNum, replace)

*re-label to make it easier to convey what's being measured
gen boundarytype = firstsyll
label variable boundarytype "boundarytype"
label define cmpsetlabel 1 "+S+M", modify
label define boundarylabel 1 "+C1+S" 0 "-C1+S"
label values boundarytype boundarylabel


codebook distsyllNum
mixed keytime i.boundarytype##i.distsyllNum || _all:R.subject || item: if keytime > 50 & keytime < 1000 & inrange(distsyllNum,1,2), mle nogroup nolrtest nofvlabel noretable

contrast i.boundarytype##i.distsyllNum
contrast i.distsyllNum@i.boundarytype

quietly margins distsyllNum#boundarytype
marginsplot,  ytitle ("Predicted latencies (ms)") xtitle("Position relative to the start of the first letter of the second syllable") title("") legend(rows(2) ring(0) pos(11)) byopts(title(" ")) subtitle(, bcolor(none) )
graph save  "output/figures/fig_E1_b_syll", replace
graph export  "output/figures/fig_E1_b_syll.png", replace
restore