04-Main-Analyses.Rmd

# Main Analyses

After having visualised the data, we will now have a look at the research questions and hypotheses. You can find a complete overview of our research questions and hypotheses in our pre-registration [https://osf.io/wuqy9/]. For all analyses, we followed the pre-registered plan for dealing with non-convergence. However, sometimes a model would convergence but afex::mixed would produce negative eigenvalues. In these cases, we attempted to remove interactions in the random slopes or the random correlations. We also checked the following assumptions of our models: 1) Outliers, 2) Homoscedasticity 3) Normality. We also attempted to check for influential cases, but the functions typically used to do that often failed to converge, even if lme4::lmer() converged. Therefore, checks for influential cases are not included in the present report.

## Research Question 1 

Research question 1 is "Do participants report lower-arousal types of affect when reporting about alone period than when reporting about social period?".


### Hypothesis 1a

 Hypothesis 1a states that there should be lower levels of high-arousal positive affect in alone compared to social periods. Therefore, we fit a mixed linear model.



```{r mixedmodel1, echo = TRUE}
lme18RQ1d <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> converges

#get first summary
summary(lme18RQ1d)

```


Next, we want to check the assumptions of our model.


```{r assumptions1, echo = TRUE}
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ1d , scaled = TRUE)) > 2) / length(resid(lme18RQ1d ))
sum(abs(resid(lme18RQ1d , scaled = TRUE)) > 2.5) / length(resid(lme18RQ1d )) 
sum(abs(resid(lme18RQ1d , scaled = TRUE)) > 3) / length(resid(lme18RQ1d )) 

#2) Homoscedasticity
plot(lme18RQ1d, type = c('p', 'smooth'))

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ1d, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ1d, scaled = TRUE))


```

After, we calculate p-values using afex::mixed.

```{r pvalue1, echo = TRUE}
RQ1d_p <- afex::mixed (HighArousalPositiveAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR", test_intercept = TRUE, cl = MyCluster)
#effect of alone significant
#positive estimate
#Alone = -1, people = 1
#get means

RQ1d_p


df_analysis18_2 <- df_analysis18 %>% select(LowArousalNegativeAffectMean, LowArousalPositiveAffectMean, HighArousalNegativeAffectMean, HighArousalPositiveAffectMean, Alone)

by(df_analysis18_2,df_analysis18_2$Alone, summary)

#check contrasts to verify results
contrasts(df_analysis18$Alone) 

```

We find support for hypothesis 1a.

### Hypothesis 1b

Hypothesis 1b was that we expect lower levels of high arousal negative affect in alone compared to social episodes.


```{r mixedmodels2, echo = TRUE}
lme18RQ1c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18) 


#get first summary
summary(lme18RQ1c)

```

Next, we again check the assumptions.

```{r assumptions2, echo = TRUE}
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ1c , scaled = TRUE)) > 2) / length(resid(lme18RQ1c )) 
sum(abs(resid(lme18RQ1c , scaled = TRUE)) > 2.5) / length(resid(lme18RQ1c )) 
sum(abs(resid(lme18RQ1c , scaled = TRUE)) > 3) / length(resid(lme18RQ1c )) 

#2) Homoscedasticity
plot(lme18RQ1c, type = c('p', 'smooth'))
#Looks mostly fine, but strays quite far at the end

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ1c, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ1c, scaled = TRUE))

#Normal enough, but the qqplot shows a lot of outliers in the tails, which is quite problematic

```

The assumptions are somewhat violated, hence interpretation of the results should be done with caution.

```{r pvalues2, echo = TRUE}
RQ1c_p <- afex::mixed (HighArousalNegativeAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR", test_intercept = TRUE, cl = MyCluster)

RQ1c_p
```

Hypothesis 1b is not confirmed.

### Hypothesis 1c

Hypothesis 1c states that there should be higher levels of low arousal positive affect in alone compared to social periods.

```{r mixedmodels3, echo = TRUE}
lme18RQ1b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18) 
#get first summary
summary(lme18RQ1b)
```

Here, we again check our assumptions.

```{r assumptions3, echo = TRUE}
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ1b , scaled = TRUE)) > 2) / length(resid(lme18RQ1b ))
sum(abs(resid(lme18RQ1b , scaled = TRUE)) > 2.5) / length(resid(lme18RQ1b )) 

sum(abs(resid(lme18RQ1b , scaled = TRUE)) > 3) / length(resid(lme18RQ1b )) 

#2) Homoscedasticity
plot(lme18RQ1b, type = c('p', 'smooth'))
#Looks mostly fine

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ1b, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ1b, scaled = TRUE))

#Normal enough, but the qqplot shows quite some outliers in the tails
```

There seem to be some minor problems with the assumptions, proceed with caution.

```{r pvalues3, echo = TRUE}
#Calculate p-values
RQ1b_p <- afex::mixed (LowArousalPositiveAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR", test_intercept = TRUE, cl = MyCluster)
#Effect of Alone is significant

RQ1b_p 

#get means
df_analysis18_2 <- df_analysis18 %>% select(LowArousalNegativeAffectMean, LowArousalPositiveAffectMean, HighArousalNegativeAffectMean, HighArousalPositiveAffectMean, Alone)
#estimate is positive, mean is lower for alone == 1, hence people have higher levels of lowarousalpositive affect when they are among people
by(df_analysis18_2,df_analysis18_2$Alone, summary)

#check contrasts to verify results
contrasts(df_analysis18$Alone) 
#Alone = -1, People = 1, 
#Positive Fixed Effect, hence it fits with the means.
#Hypothesis 1C not confirmed, effect in the opposite direction.
```

Hypothesis 1c is not confirmed. In contrast, the effect seems to go in the opposite direction.

### Hypothesis 1d

Hypothesis 1d states that there should be higher levels of low-arousal negative affect in alone compared to social periods.

```{r mixedmodels4, echo = TRUE}
lme18RQ1 <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18)

#get first summary
summary(lme18RQ1)




#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ1 , scaled = TRUE)) > 2) / length(resid(lme18RQ1 )) 
sum(abs(resid(lme18RQ1 , scaled = TRUE)) > 2.5) / length(resid(lme18RQ1 )) 
sum(abs(resid(lme18RQ1 , scaled = TRUE)) > 3) / length(resid(lme18RQ1 )) 

#2) Homoscedasticity
plot(lme18RQ1, type = c('p', 'smooth'))
#Also somewhat violated

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ1, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ1, scaled = TRUE))

#Normal enough, but the qqplot shows quite some outliers in the tails
```


The assumptions seem somewhat violated.

```{r pvalues4, echo = TRUE}
RQ1_p <- afex::mixed (LowArousalNegativeAffectMean ~ Alone + (1 + Alone | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR", test_intercept = TRUE, cl = MyCluster)

RQ1_p



#Hypothesis 1D is supported.

by(df_analysis18_2,df_analysis18_2$Alone, summary)

#check contrasts to verify results
contrasts(df_analysis18$Alone) 
```

Hypothesis 1d is supported.

### Correction for Multiple Testing

Now, we need to check whether 1a and 1d remain significant even after correction for multiple testing. We pre-registered the Holm correction.

```{r mtrq1, echo = TRUE}
#Check if analyses still significant after correction

df_analysis18RQ1_pvalues <- as.data.frame(rbind(RQ1_p$anova_table$`Pr(>F)`, RQ1b_p$anova_table$`Pr(>F)`, RQ1c_p$anova_table$`Pr(>F)`, RQ1d_p$anova_table$`Pr(>F)`))

df_analysis18RQ1_pvalues$adjusted <- p.adjust(df_analysis18RQ1_pvalues$V2, method = "holm")

df_analysis18RQ1_pvalues$adjusted 
```

The p-value associated with Hypothesis 1a and 1d remains significant. Moreover, the p-value associated with the effect opposite to the predicted effect in 1c remains significant.

## Research Question 2

Research Question 2 is: "Does day-level exposure to nature correlate with affect in that instance, or moderate the above effect?".

### Hypotheses 2Aa and 2bc

2Aa: We expect a positive correlation between day-level exposure to nature and high arousal positive affect. 

2Bc: We have no prediction that Day-level exposure to nature will moderate the difference between Alone vs. Social on high-arousal types of affect.

```{r 2Aa/2bc, echo = TRUE }
lme18RQ2d <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ2d <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ2d <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges



#get first summary
summary(lme18RQ2d)
```

Next, we check assumptions:

```{r 2Aa assumptions, echo = TRUE}
#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ2d , scaled = TRUE)) > 2) / length(resid(lme18RQ2d ))
sum(abs(resid(lme18RQ2d , scaled = TRUE)) > 2.5) / length(resid(lme18RQ2d )) 
sum(abs(resid(lme18RQ2d , scaled = TRUE)) > 3) / length(resid(lme18RQ2d )) 
#2) Homoscedasticity
plot(lme18RQ2d, type = c('p', 'smooth'))
#looks mostly fine

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ2d, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ2d, scaled = TRUE))
#qqplot shows quite some outliers in the tails

```

As usual, there are some issues with the assumptions here.

```{r 2Aap, echo = TRUE}
RQ2d_p <- afex::mixed (HighArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)

RQ2d_p

#2Aa not supported. "Support" for 2Bc/d because we predicted no effect and there is no effect
#Significant effect of Alone
```

There is no support for Hypothesis 2Aa. However, as predicted, there is no moderating effect (i.e., 2Bc confirmed).

### Hypotheses 2Ab and 2Bd

2Ab: We expect that day-level exposure to nature will negatively correlate with high-arousal negative affect.

2Bd: We have no prediction that Day-level exposure to nature will moderate the difference between Alone vs. Social on high-arousal types of affect.


```{r 2Ab, echo = TRUE}

lme18RQ2c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ2c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ2c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges 


#get first summary
summary(lme18RQ2c)
```

Next, we check the assumptions

```{r 2Abassumptions, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ2c , scaled = TRUE)) > 2) / length(resid(lme18RQ2c ))
sum(abs(resid(lme18RQ2c , scaled = TRUE)) > 2.5) / length(resid(lme18RQ2c ))
sum(abs(resid(lme18RQ2c , scaled = TRUE)) > 3) / length(resid(lme18RQ2c ))

#2) Homoscedasticity
plot(lme18RQ2c, type = c('p', 'smooth'))
#looks mostly fine

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ2c, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ2c, scaled = TRUE))

#Barely Normal and the qqplot shows quite some outliers in the tails
```

Again, there are quite some problems with the assumptions. Proceed with caution.

```{r 2Abp, echo = TRUE}
RQ2c_p <- afex::mixed (HighArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)

RQ2c_p

#2Ab not supported. 2Bc/d "supported", as in we did not predict an effect and we did not observe an effect
```

We found no support for hypothesis 2Ab, but we did find the predicted absence of a moderating effect (2Bd)

### Hypotheses 2Ac and 2Ba

2Ac: We expect day-level exposure to nature to corelative positive with low-arousal positive affect.
2Ba: Day-level exposure to nature will moderate 
the difference between Alone vs. Social on low-arousal positive affect.

```{r 2Ac, echo = TRUE}
lme18RQ2b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ2b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ2b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) 


#get first summary
summary(lme18RQ2b)
```

Next, we again check the assumptions.

```{r 2Ac assumptions, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ2b , scaled = TRUE)) > 2) / length(resid(lme18RQ2b )) 
sum(abs(resid(lme18RQ2b , scaled = TRUE)) > 2.5) / length(resid(lme18RQ2b )) 
sum(abs(resid(lme18RQ2b , scaled = TRUE)) > 3) / length(resid(lme18RQ2b )) 

#2) Homoscedasticity
plot(lme18RQ2b, type = c('p', 'smooth'))
#looks mostly fine

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ2b, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ2b, scaled = TRUE))

#Normal enough, but the qqplot shows quite some outliers in the tails
```

Minor problems, but better than the preceeding models.

```{r 2Acp, echo = TRUE}
RQ2b_p <- afex::mixed (LowArousalPositiveAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)

RQ2b_p 

#2Ac and 2Ba not supported
```

We find no support for hypothesis 2Ac or hypothesis 2Ba.


### Hypotheses 2Ad and 2Bb
2Ad: We expect day-level exposure to nature to negatively correlate with low-arousal negative affect.

2Bb: Day-level exposure to nature will moderate the difference between Alone vs. Social on low-arousal negative affect.

```{r 2Ad, echo = TRUE}

lme18RQ2 <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ2 <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ2 <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges

#get first summary
summary(lme18RQ2)
```

Next, we check the assumptions.

```{r 2Ad ass, echo = TRUE}
#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ2 , scaled = TRUE)) > 2) / length(resid(lme18RQ2 )) 
sum(abs(resid(lme18RQ2 , scaled = TRUE)) > 2.5) / length(resid(lme18RQ2 )) 
sum(abs(resid(lme18RQ2 , scaled = TRUE)) > 3) / length(resid(lme18RQ2 ))

#2) Homoscedasticity
plot(lme18RQ2, type = c('p', 'smooth'))
#Trails off strongly in the tails

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ2, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ2, scaled = TRUE))

#Normal enough, but the qqplot shows quite some outliers in the tails
```

As usual, some problems with the assumptions.

```{r 2Adp, echo = TRUE}
RQ2_p <- afex::mixed (LowArousalNegativeAffectMean ~ Alone*DayLevelNatureScaled + (1 + Alone*DayLevelNatureScaled | Participant.ID) + (1 | SubjDay), control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), data = df_analysis18, type = 3, method = "KR", test_intercept = TRUE, cl = MyCluster)

RQ2_p 

#no significant hypothesized effects, 
#Hypothesis 2Ad or 2Bb not supported
```

We find no support for hypothesis 2Ad and hypothesis 2Bb.

## Research Question 3
Research Question 3 is "Does day-level intentionality for solitude correlate with affect in that instance, or moderate the above effect?".

### Hypotheses 3Aa and 3Bc
3Aa: We expect a positive correlation between day-level intentionality for solitude and high arousal positive affect.

3Bc: We will explore the interaction but 
have no prediction that day-level solitude intentionality will moderate the difference between Alone vs. Social on high-arousal types of affect.

```{r 3Aa, echo = TRUE}
lme18RQ3d <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #converges

lme18RQ3d <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa"))



#get first summary
summary(lme18RQ3d)
```

Next, we check the assumptions.

```{r 3Aaassum, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ3d , scaled = TRUE)) > 2) / length(resid(lme18RQ3d )) 
sum(abs(resid(lme18RQ3d , scaled = TRUE)) > 2.5) / length(resid(lme18RQ3d )) 
sum(abs(resid(lme18RQ3d , scaled = TRUE)) > 3) / length(resid(lme18RQ3d ))

#2) Homoscedasticity
plot(lme18RQ3d, type = c('p', 'smooth'))
#looks fine
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ3d, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ3d, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

We again see some issues with the assumptions here.

```{r 3Aap, echo = TRUE}
RQ3d_p <- afex::mixed (HighArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster, control = lmerControl(optimizer = "bobyqa"))
#3Aa not supported. 3Bc/d: No interaction
RQ3d_p 
```

We find no support for hypothesis 3Aa. However, we do find the predicted absence of a moderating effect (3Bc).

### Hypotheses 3Ab and 3Bd
3Ab: We expect that day-level intentionality for solitude will negatively correlate with high arousal negative affect.

3Bd: We will explore the interaction but 
have no prediction that day-level solitude intentionality will moderate the difference between Alone vs. Social on high-arousal types of affect.

```{r 3Ab, echo = TRUE}

lme18RQ3c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ3c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ3c <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges




#get first summary
summary(lme18RQ3c)
```

Now, we again test the assumptions.

```{r 3Abassum, echo = TRUE}
#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ3c , scaled = TRUE)) > 2) / length(resid(lme18RQ3c ))
sum(abs(resid(lme18RQ3c , scaled = TRUE)) > 2.5) / length(resid(lme18RQ3c ))
sum(abs(resid(lme18RQ3c , scaled = TRUE)) > 3) / length(resid(lme18RQ3c )) 

#2) Homoscedasticity
plot(lme18RQ3c, type = c('p', 'smooth'))
#looks violated

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ3c, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ3c, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

As usual, we see that there are some issues with the assumptions.


```{r 3Abp, echo = TRUE}
RQ3c_p <- afex::mixed (HighArousalNegativeAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#3Ab not supported, 3Bc/d: No interaction

RQ3c_p 
```

We find no support for hypothesis 3Ab. However, we do find the predicted absence of a moderating effect (3Bd)

### Hypotheses 3Ac and 3Ba.

3Ac: We expect a positive correlation between day-level intentionality for solitude and low-arousal positive affect.

3Ba: Day-level solitude intentionality will moderate the difference between Alone vs. Social on low-arousal positive affect.


```{r 3Ac, echo = TRUE}
lme18RQ3b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ3b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ3b <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE))


#get first summary
summary(lme18RQ3b)
```

Next, we check our assumptions.

```{r 3Acassum, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ3b , scaled = TRUE)) > 2) / length(resid(lme18RQ3b ))
sum(abs(resid(lme18RQ3b , scaled = TRUE)) > 2.5) / length(resid(lme18RQ3b )) 
sum(abs(resid(lme18RQ3b , scaled = TRUE)) > 3) / length(resid(lme18RQ3b ))

#2) Homoscedasticity
plot(lme18RQ3b, type = c('p', 'smooth'))
#looks good

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ3b, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ3b, scaled = TRUE))

#qqplot shows some outliers in the tails
```

As usual, assumptions are somewhat violated.

```{r 3Acp, echo = TRUE}
RQ3b_p <- afex::mixed (LowArousalPositiveAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#3Ac and 3Ba not supported

RQ3b_p
```

We find no support for hypothesis 3Ac or hypothesis 3Ba.

### Hypotheses 3Ad and 3Bb

3Ad: We expect that day-level intentionality for solitude negatively correlates with low arousal negative affect.

3Bb: Day-level solitude intentionality will moderate the difference between Alone vs. Social 
on low-arousal negative affect.

```{r 3Ad, echo = TRUE}
lme18RQ3 <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) 

#get first summary
summary(lme18RQ3)

```

Next, we test our assumptions.

```{r 3Adassum, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ3b , scaled = TRUE)) > 2) / length(resid(lme18RQ3b ))
sum(abs(resid(lme18RQ3b , scaled = TRUE)) > 2.5) / length(resid(lme18RQ3b )) 
sum(abs(resid(lme18RQ3b , scaled = TRUE)) > 3) / length(resid(lme18RQ3b ))
#2) Homoscedasticity
plot(lme18RQ3b, type = c('p', 'smooth'))
#looks good

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ3b, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ3b, scaled = TRUE))

#qqplot shows some outliers in the tails
```

Assumptions are slightly violated.

```{r 3Adp, echo = TRUE}
RQ3b_p <- afex::mixed (LowArousalNegativeAffectMean ~ Alone*DayLevelIntentionalitySolitudeScaled + (1 + Alone*DayLevelIntentionalitySolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#3Ac and 3Ba not supported
RQ3b_p

```

## Research Question 4

Research question 4 is: "Does structure for solitude (at both individual levels and day-level) correlate with affect at the event level or moderate the above effect?"

First, we will run the models investigating day-level associations.

### Hypotheses 4Aa and 4Ba Day-Level

4Aa: We expect day-level structure for solitude to positively correlate with high-arousal positive affect.
4Ba: Day-level/Person-level structure for solitude will moderate the difference between Alone vs. Social on high-arousal positive affect.

```{r 4Aa day, echo = TRUE}
lme18RQ4d.day <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) 

#get first summary
summary(lme18RQ4d.day)
```

Next, we check the assumptions.

```{r 4Aaday assumptions, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4d.day , scaled = TRUE)) > 2) / length(resid(lme18RQ4d.day )) 
sum(abs(resid(lme18RQ4d.day , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4d.day )) 
sum(abs(resid(lme18RQ4d.day , scaled = TRUE)) > 3) / length(resid(lme18RQ4d.day ))

#2) Homoscedasticity
plot(lme18RQ4d.day, type = c('p', 'smooth'))
#looks good

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4d.day, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4d.day, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

Again, the assumptions are slightly violated.

```{r 4Aadayp, echo = TRUE}
RQ4d.day_p <- afex::mixed (HighArousalPositiveAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#Hypothesis 4Aa and 4Ba day not supported
RQ4d.day_p
```

We find no support for hypothesis 4Aa and 4Ba on a day-level.

### Hypotheses 4Ab and 4Bb Day-level

4Ab: We expect day-level structure for solitude to correlate negatively with high arousal negative affect.

4Bb: Day-level/Person-level structure for solitude 
will moderate the difference between Alone vs. Social on high-arousal negative affect.

```{r 4Ab day, echo = TRUE}

lme18RQ4c.day <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ4c.day <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #does not convergence

lme18RQ4c.day <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges



#get first summary
summary(lme18RQ4c.day)
```

Let´s check the assumptions again:

```{r 4Ab day assum, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4c.day , scaled = TRUE)) > 2) / length(resid(lme18RQ4c.day )) 
sum(abs(resid(lme18RQ4c.day , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4c.day )) 
sum(abs(resid(lme18RQ4c.day , scaled = TRUE)) > 3) / length(resid(lme18RQ4c.day )) 

#2) Homoscedasticity
plot(lme18RQ4c.day, type = c('p', 'smooth'))
#looks good

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4c.day, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4c.day, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

The assumptions are again slightly violated.

```{r 4Ab day p, echo = TRUE}
RQ4c.day_p <- afex::mixed (HighArousalNegativeAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#Hypothesis 4Ab and 4Bb day not supported
RQ4c.day_p 
```

Hypothesis 4Ab and 4Bb on a day-level are not supported.

### Hypotheses 4Ac and 4Bc Day Level

4Ac: We expect day-level structure for solitude to positively correlate with low-arousal positive affect.

4Bc: We will look at the interaction but 
have no prediction that day-level/person-level structure for solitude willmoderate the difference between Alone vs. Social on low-arousal types of affect.

```{r 4Ac day, echo = TRUE}
lme18RQ4b.day <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> not converges

lme18RQ4b.day <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18,  control = lmerControl(optimizer = "bobyqa")) #-> not

lme18RQ4b.day <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) -> converges


#get first summary
summary(lme18RQ4b.day)
```

Let´s check the assumptions next.

```{r 4Ac day assum, echo = TRUE}
#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4b.day , scaled = TRUE)) > 2) / length(resid(lme18RQ4b.day ))
sum(abs(resid(lme18RQ4b.day , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4b.day ))
sum(abs(resid(lme18RQ4b.day , scaled = TRUE)) > 3) / length(resid(lme18RQ4b.day ))

#2) Homoscedasticity
plot(lme18RQ4b.day, type = c('p', 'smooth'))
#looks fine

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4b.day, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4b.day, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

As usual, the assumptions are somewhat violated.

```{r 4Ac day p, echo = TRUE}
RQ4b.day_p <- afex::mixed (LowArousalPositiveAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#Hypothesis 4Ac day not supported, 4Bc/d: No Interaction
RQ4b.day_p

```

We find no evidence for hypothesis 4Ac on a day-level. Moreover, as predicted, we find no moderating effect (4Bc).

### Hypotheses 4Ad and 4Bd Day-Level

4Ad: We expect day-level structure for solitude to negatively correlate with low-arousal negative affect.

4Bd: We will look at the interaction but 
have no prediction that day-level/person-level structure for solitude willmoderate the difference between Alone vs. Social on low-arousal types of affect.

```{r 4Ad day}
lme18RQ4.day <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) 


#get first summary
summary(lme18RQ4.day)

```

Next, we again test our assumptions.

```{r 4Ad day assum, echo = TRUE}
#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4.day , scaled = TRUE)) > 2) / length(resid(lme18RQ4.day ))
sum(abs(resid(lme18RQ4.day , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4.day ))
sum(abs(resid(lme18RQ4.day , scaled = TRUE)) > 3) / length(resid(lme18RQ4.day ))

#2) Homoscedasticity
plot(lme18RQ4.day, type = c('p', 'smooth'))
#looks violated

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4.day, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4.day, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

As usual, assumptions are somewhat violated.

```{r 4Ad day p, echo = TRUE}
RQ4.day_p <- afex::mixed (LowArousalNegativeAffectMean ~ Alone*DayLevelStructureSolitudeScaled + (1 + Alone*DayLevelStructureSolitudeScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
RQ4.day_p 
#Hypothesis 4Ad day not supported, 4Bc/d: No Interaction
```

We find no evidence for hypothesis 4Ad. Moreover, as expected, we find no moderating effect (4Bd).

### Hypothesis 4Aa and 4Ba Person-Level

4Aa: We expect person-level structure for solitude to psotiviely correlate with high-arousal positive affect.

4Ba: Day-level/Person-level structure for solitude will moderate the difference between Alone vs. Social on high-arousal positive affect

```{r 4Aa person, echo = TRUE}
lme18RQ4d.individual <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ4d.individual <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ4d.individual <- lme4::lmer(HighArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE))

#get first summary
summary(lme18RQ4d.individual)

```

Next, we check our assumptions.

```{r 4Aa person assum, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4d.individual , scaled = TRUE)) > 2) / length(resid(lme18RQ4d.individual ))
sum(abs(resid(lme18RQ4d.individual , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4d.individual )) 
sum(abs(resid(lme18RQ4d.individual , scaled = TRUE)) > 3) / length(resid(lme18RQ4d.individual ))

#2) Homoscedasticity
plot(lme18RQ4d.individual, type = c('p', 'smooth'))
#looks good
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4d.individual, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4d.individual, scaled = TRUE))

#qqplot shows outliers in the tails
```

As usual, there is some violation of the assumptions.

```{r 4Aa person p, echo = TRUE}
#RQ4d.individual_p <- afex::mixed (HighArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)

#Fehler in h(simpleError(msg, call)) : Fehler bei der Auswertung des Argumentes 'x' bei der Methodenauswahl für Funktion 'forceSymmetric': System ist für den Rechner singulär: reziproke Konditionszahl = 3.16421e-35

RQ4d.individual_p <- afex::mixed (HighArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "S", cl = MyCluster)

RQ4d.individual_p
#Hypothesis 4Aa and 4Ba not suported

```

Hypothesis 4Aa is not supported. Hypothesis 4Ba is supported pending correction for multiple testing.

### Hypotheses 4Ab and 4Bb Person-Level

4Ab: We expect person-level structure for solitude to be negatively correlated with high-arousal negative affect.

4Bb: Day-level/Person-level structure for solitude 
will moderate the difference between Alone vs. 
Social on high-arousal negative affect.

```{r 4Ab person}
lme18RQ4c.individual <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ4c.individual <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ4c.individual <- lme4::lmer(HighArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) 

#get first summary
summary(lme18RQ4c.individual)

```

Next, we check the assumptions.

```{r 4Ab person assump, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4c.individual , scaled = TRUE)) > 2) / length(resid(lme18RQ4c.individual ))
sum(abs(resid(lme18RQ4c.individual , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4c.individual ))
sum(abs(resid(lme18RQ4c.individual , scaled = TRUE)) > 3) / length(resid(lme18RQ4c.individual ))

#2) Homoscedasticity
plot(lme18RQ4c.individual, type = c('p', 'smooth'))
#looks slightly violated

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4c.individual, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4c.individual, scaled = TRUE))

#qqplot shows strong outliers in the tails
```

There seem to be quite some violations there. 

```{r 4Ab person p, echo = TRUE}
#RQ4c.individual_p <- afex::mixed (HighArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)

#Fehler in h(simpleError(msg, call)) : Fehler bei der Auswertung des Argumentes 'x' bei der Methodenauswahl für Funktion 'forceSymmetric': Lapackroutine dgesv: System ist genau singulär: U[12,12] = 0


RQ4c.individual_p <- afex::mixed (HighArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "S",  cl = MyCluster)

#Hypothesis 4Ab individual not supported, 4Bb There is an interaction
RQ4c.individual_p 
```

Hypothesis 4Ab and hypothesis 4Bb are not supported.


### Hypotheses 4Ac and 4Bc Person-Level

4Ac: We expect person-level structure for solitude to positively correlate with low-arousal positive affect.

4Bc: We will look at the interaction but have no prediction that day-level/person-level structure for solitude will moderate the difference between Alone vs. Social on low-arousal types of affect.

```{r 4Ac person, echo = TRUE}
lme18RQ4b.individual <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ4b.individual <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ4b.individual <- lme4::lmer(LowArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges



#get first summary
summary(lme18RQ4b.individual)

```

Next, we test our assumptions.

```{r 4Ac person assump, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4b.individual , scaled = TRUE)) > 2) / length(resid(lme18RQ4b.individual )) 
sum(abs(resid(lme18RQ4b.individual , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4b.individual ))
sum(abs(resid(lme18RQ4b.individual , scaled = TRUE)) > 3) / length(resid(lme18RQ4b.individual )) 

#2) Homoscedasticity
plot(lme18RQ4b.individual, type = c('p', 'smooth'))
#looks fine
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4b.individual, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4b.individual, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

As usual, assumptions are somewhat violated.

```{r 4Ac person p, echo = TRUE}
#RQ4b.individual_p <- afex::mixed (LowArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
#produces error 

RQ4b.individual_p <- afex::mixed (LowArousalPositiveAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "S",  cl = MyCluster)
#Hypothesis 4Ac individual not supported, 4Bc/d: No Interaction
RQ4b.individual_p 
```

Hypothesis 4Ac is not supported. Moreover, contrary to our predictions there is a moderating effect.

### Hypotheses 4Ad and 4Bd Person-Level

4Ad: We expect person-level structure for solitude to negatively correlate with low-arousal negative affect.

4Bd: We will look at the interaction but have no prediction that day-level/person-level structure for solitude will moderate the difference between Alone vs. Social on low-arousal types of affect.

```{r 4Ad person, echo = TRUE}
lme18RQ4.individual <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ4.individual <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ4.individual <- lme4::lmer(LowArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE)) #-> converges
#get first summary
summary(lme18RQ4.individual)
```

Next, we test our assumptions.

```{r 4Ad person assump, echo = TRUE}

#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ4.individual , scaled = TRUE)) > 2) / length(resid(lme18RQ4.individual )) 
sum(abs(resid(lme18RQ4.individual , scaled = TRUE)) > 2.5) / length(resid(lme18RQ4.individual )) 
sum(abs(resid(lme18RQ4.individual , scaled = TRUE)) > 3) / length(resid(lme18RQ4.individual )) 

#2) Homoscedasticity
plot(lme18RQ4.individual, type = c('p', 'smooth'))
#looks violated

#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ4.individual, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ4.individual, scaled = TRUE))

#qqplot shows quite some outliers in the tails
```

As usual, there are some issues with the assumptions here.

```{r 4Ad person p, echo = TRUE}
#RQ4.individual_p <- afex::mixed (LowArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)

#Fehler in h(simpleError(msg, call)) : Fehler bei der Auswertung des Argumentes 'x' bei der Methodenauswahl für Funktion 'forceSymmetric': System ist für den Rechner singulär: reziproke Konditionszahl = 1.36756e-33

RQ4.individual_p <- afex::mixed (LowArousalNegativeAffectMean ~ Alone*StructureForSolitudeIndividualScaled + (1 + Alone*StructureForSolitudeIndividualScaled  | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE), type = 3, method = "S", cl = MyCluster)


#Hypothesis 4Ad individual not supported, 4Bc/d: No Interaction
RQ4.individual_p 
```

Hypothesis 4Ad is not supported. However, as predicted, there is no moderating effect.


### Correction for Multiple Testing

Now, we need to check whether 1a remains significant even after correction for multiple testing. We pre-registered the Holm correction.

```{r mtrq4, echo = TRUE}
#Check if analyses still significant after correction

df_analysis18RQ4_person_pvalues <- as.data.frame(rbind(RQ4.individual_p$anova_table$`Pr(>F)`[3], RQ4b.individual_p$anova_table$`Pr(>F)`[3], RQ4c.individual_p$anova_table$`Pr(>F)`[3], RQ4d.individual_p$anova_table$`Pr(>F)`[3]))

df_analysis18RQ4_person_pvalues$adjusted <- p.adjust(df_analysis18RQ4_person_pvalues$V1, method = "holm")

df_analysis18RQ4_person_pvalues$adjusted 
```

After correction for multiple testing all significant results disappear.


## Exploratory Research Questions

### Preparation

Finally, we perform a pre-registered exploratory analysis without pre-defined hypothesis. This corresponds to RQ5 in the pre-registration: " Is affective levels of an alone period linked to  affective levels of an adjacent social period, when people switch from one experience to another, depending on which one comes first?"

First, we create the additional variables we require for this analysis.

```{r RQ5prep, echo = TRUE}

#create lagged alone variable as preparation for creating the epoch variable
df_analysis18 <- 
  df_analysis18 %>%
  group_by(Participant.ID) %>%
  mutate(AloneLagged = dplyr::lag(Alone, n = 1, default = NA))

#create epoch. epoch = 1 if alone & alone lagged are not equal and not NA.
df_analysis18 <- 
  df_analysis18 %>%
  mutate(epoch = ifelse(Alone == 1 & AloneLagged == 0 | Alone == 0 & AloneLagged == 1, 1, 0)) 

df_analysis18$epoch <- as.factor(as.character(df_analysis18$epoch))

#create affect.0 variables, which are the lagged versions of the affect variables.
df_analysis18 <- 
  df_analysis18 %>%
  group_by(Participant.ID) %>%
  mutate(LowArousalNegativeAffectMeanLagged = dplyr::lag(LowArousalNegativeAffectMean, n = 1, default = NA))

df_analysis18 <- 
  df_analysis18 %>%
  group_by(Participant.ID) %>%
  mutate(LowArousalPositiveAffectMeanLagged = dplyr::lag(LowArousalPositiveAffectMean, n = 1, default = NA))

df_analysis18 <- 
  df_analysis18 %>%
  group_by(Participant.ID) %>%
  mutate(HighArousalNegativeAffectMeanLagged = dplyr::lag(HighArousalNegativeAffectMean, n = 1, default = NA))

df_analysis18 <- 
  df_analysis18 %>%
  group_by(Participant.ID) %>%
  mutate(HighArousalPositiveAffectMeanLagged = dplyr::lag(HighArousalPositiveAffectMean, n = 1, default = NA))

#scale the lagged variable as predictor

df_analysis18$LowArousalNegativeAffectMeanLaggedScaled <- scale(df_analysis18$LowArousalNegativeAffectMeanLagged)

df_analysis18$LowArousalPositiveAffectMeanLaggedScaled <- scale(df_analysis18$LowArousalPositiveAffectMeanLagged)

df_analysis18$HighArousalNegativeAffectMeanLaggedScaled <- scale(df_analysis18$HighArousalNegativeAffectMeanLagged)

df_analysis18$HighArousalPositiveAffectMeanLaggedScaled <- scale(df_analysis18$HighArousalPositiveAffectMeanLagged)
```


### Analyses

Next, we perform all four analysis and perform a correction for multiple testing where necessary.

```{r RQ5analysis, echo = TRUE}

#####fit mixed model#####
##LowArousalNegative Affect##
lme18RQ5 <- lme4::lmer(LowArousalNegativeAffectMean ~ LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence


lme18RQ5 <- lme4::lmer(LowArousalNegativeAffectMean ~ LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) # does not converge

lme18RQ5 <- lme4::lmer(LowArousalNegativeAffectMean ~ LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5))) 

#get first summary
summary(lme18RQ5)



#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ5 , scaled = TRUE)) > 2) / length(resid(lme18RQ5 ))
sum(abs(resid(lme18RQ5 , scaled = TRUE)) > 2.5) / length(resid(lme18RQ5 )) 
sum(abs(resid(lme18RQ5 , scaled = TRUE)) > 3) / length(resid(lme18RQ5 )) 

#2) Homoscedasticity
plot(lme18RQ5, type = c('p', 'smooth'))
#looks good
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ5, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ5, scaled = TRUE))

#qqplot shows outliers in the tails


#Conclusion: There are some issues with the assumptions here. Proceed with caution.

#Calculate p-values
#says the predictor is not centered on zero
RQ5_p <- afex::mixed (LowArousalNegativeAffectMean ~ LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalNegativeAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5)), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
RQ5_p

##low arousal positive affect##
#####fit mixed model#####
lme18RQ5b <- lme4::lmer(LowArousalPositiveAffectMean ~ LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ5b <- lme4::lmer(LowArousalPositiveAffectMean ~ LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ5b <- lme4::lmer(LowArousalPositiveAffectMean ~ LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5))) 

#get first summary
summary(lme18RQ5b)



#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ5b , scaled = TRUE)) > 2) / length(resid(lme18RQ5b ))
sum(abs(resid(lme18RQ5b , scaled = TRUE)) > 2.5) / length(resid(lme18RQ5b )) 
sum(abs(resid(lme18RQ5b , scaled = TRUE)) > 3) / length(resid(lme18RQ5b )) 

#2) Homoscedasticity
plot(lme18RQ5b, type = c('p', 'smooth'))
#looks good
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ5b, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ5b, scaled = TRUE))

#qqplot shows outliers in the tails


#Conclusion: There are some issues with the assumptions here. Proceed with caution.

#Calculate p-values

RQ5b_p <- afex::mixed (LowArousalPositiveAffectMean ~ LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5)), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
RQ5b_p


##high arousal negative affect##
#####fit mixed model#####
lme18RQ5c <- lme4::lmer(HighArousalNegativeAffectMean ~ HighArousalNegativeAffectMeanLaggedScaled*Alone*epoch + (1 +HighArousalNegativeAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ5c <- lme4::lmer(HighArousalNegativeAffectMean ~ HighArousalNegativeAffectMeanLaggedScaled*Alone*epoch + (1 +HighArousalNegativeAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ5c <- lme4::lmer(LowArousalPositiveAffectMean ~ LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5))) 

#get first summary
summary(lme18RQ5c)



#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ5c , scaled = TRUE)) > 2) / length(resid(lme18RQ5c ))
sum(abs(resid(lme18RQ5c , scaled = TRUE)) > 2.5) / length(resid(lme18RQ5c )) 
sum(abs(resid(lme18RQ5c , scaled = TRUE)) > 3) / length(resid(lme18RQ5c )) 

#2) Homoscedasticity
plot(lme18RQ5c, type = c('p', 'smooth'))
#looks good
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ5c, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ5c, scaled = TRUE))

#qqplot shows outliers in the tails


#Conclusion: There are some issues with the assumptions here. Proceed with caution.

#Calculate p-values

RQ5c_p <- afex::mixed(LowArousalPositiveAffectMean ~ LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +LowArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5)), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
RQ5c_p


##high arousal positive affect##

#####fit mixed model#####
lme18RQ5d <- lme4::lmer(HighArousalPositiveAffectMean ~ HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18) #-> does not convergence

lme18RQ5d <- lme4::lmer(HighArousalPositiveAffectMean ~ HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa")) #-> does not convergence

lme18RQ5d <- lme4::lmer(HighArousalPositiveAffectMean ~ HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5))) 

#get first summary
summary(lme18RQ5d)



#check assumptions
#1) Outliers
#proportion of residuals larger than +/- 2, 2.5, 3.
sum(abs(resid(lme18RQ5d , scaled = TRUE)) > 2) / length(resid(lme18RQ5d ))
sum(abs(resid(lme18RQ5d , scaled = TRUE)) > 2.5) / length(resid(lme18RQ5d )) 
sum(abs(resid(lme18RQ5d , scaled = TRUE)) > 3) / length(resid(lme18RQ5d )) 

#2) Homoscedasticity
plot(lme18RQ5d, type = c('p', 'smooth'))
#looks good
#3) Normality
#densityplot of the scaled residuals
densityplot(resid(lme18RQ5d, scaled = TRUE))
#qqplot of the scaled residuals
car::qqPlot(resid(lme18RQ5d, scaled = TRUE))

#qqplot shows outliers in the tails


#Conclusion: There are some issues with the assumptions here. Proceed with caution.

#Calculate p-values

RQ5d_p <- afex::mixed(HighArousalPositiveAffectMean ~ HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch + (1 +HighArousalPositiveAffectMeanLaggedScaled*Alone*epoch | Participant.ID) + (1 | SubjDay), data = df_analysis18, control = lmerControl(optimizer = "bobyqa", calc.derivs = FALSE, optCtrl=list(maxfun=1e5)), type = 3, method = "KR",test_intercept = TRUE,  cl = MyCluster)
RQ5d_p
```

For all four analysis, Alone and the lagged affect variable emerge as significant predictors. Here, we are mainly interested in the effect of the lagged affect variables. Therefore, we correct for multiple testing for the effect of the lagged affect variables

### Correction for Multiple Testing

```{r RQ5mt, echo = TRUE}

#Adjust effect of affect.0 variables
df_analysis18RQ5.2_pvalues <- as.data.frame(rbind(RQ5_p$anova_table$`Pr(>F)`[2], RQ5b_p$anova_table$`Pr(>F)`[2], RQ5c_p$anova_table$`Pr(>F)`[2], RQ5d_p$anova_table$`Pr(>F)`[2]))

df_analysis18RQ5.2_pvalues$adjusted <- p.adjust(df_analysis18RQ5.2_pvalues$V1, method = "holm")
df_analysis18RQ5.2_pvalues$adjusted #all remain significant
```

All four p-values remain significant even after correction for multiple testing.