library(nlme)
library(dyadr)
riggsp <- read.csv("riggsp.csv") A run examining the effect of own Attachment Anxiety (Anxiety_A) or actor and the effect of partner’s Attachment Anxiety (Anxiety_P) on satisfaction (Sat_A), allowing for nonindependence by correlating the errors of the two members and treating dyad members as indistinguishable:
apim = gls(Sat_A ~ Anxiety_A + Anxiety_P,
na.action=na.omit,
correlation=corCompSymm (form=~1|Dyad),
data=riggsp)
smallsummary(apim)## Correlation structure of class corCompSymm representing
## Rho
## 0.6166651
##
## Residual standard error: 6.71144
##
## Value Std.Error t-value p-value
## (Intercept) 53.3388 1.4945 35.6903 0
## Anxiety_A -1.7891 0.2944 -6.0769 0
## Anxiety_P -1.2770 0.2944 -4.3376 0## 2.5 % 97.5 %
## (Intercept) 50.4097 56.268
## Anxiety_A -2.3661 -1.212
## Anxiety_P -1.8540 -0.700This is the same code as yesterday, but the key is that the same variable in included “twice” as a predictor, once for the respondent (Anxiety_A) and once for the partner (Anxiety_P).
First look at the coefficients:
Intercept = 53.33880: The predicted level of Satisfaction for those who score zero on all the predictors. Because Anxiety is not centered this value is not meaningful.
The Actor Effect of Anxiety = -1.789: The effect of a one unit change of Anxiety_A on Sat_A holding Anxiety_P constant. This result is statistically significant.
The Partner Effect of Anxiety = -1.277: The effect of a one unit change of Anxiety_P on Sat_A holding Anxiety_A constant. This result is statistically significant.
Note that actor and partner effects are nearly equal with the partner effect a bit smaller than the actor effect. More on that when patterns are discussed.
For the 95% confidence intervals:
Note that zero is not included in the 95 percent confidence interval.
To compute effect sizes, we use partial correlations because the X variables are interval.
library(ppcor)
print("Anxiety actor effect effect size:")## [1] "Anxiety actor effect effect size:" as.numeric(pcor.test(riggsp$Anxiety_A,riggsp$Sat_A,riggsp$Anxiety_P)[1])## [1] -0.3112176print ("Anxiety partner effect effect size:") ## [1] "Anxiety partner effect effect size:" as.numeric(pcor.test(riggsp$Anxiety_P,riggsp$Sat_A,riggsp$Anxiety_A)[1])## [1] -0.2276168The effect size for actor is -.311, a moderate effect, and for partner -.228, a small effect.
The standard deviation of the errors, the unexplained variation, is 6.71144.
The multiple correlation squared can be computed by comparing the error variance for the model to a model with no predictors, which is called the Empty model. (Sometimes this is called a pseudo R squared.)
# Empty Model
apimie = summary(gls(Sat_A ~ 1,
na.action=na.omit,
correlation=corCompSymm (form=~1|Dyad),
data=riggsp))
# sd of errors for the model or esd
esd = as.numeric(apim[6])
# sd of errors for the empty model or esd0
esd0 = as.numeric(apimie[6])
# the R squared, using the crsp function
crsp (esd,esd0)## [1] 0.1582226The function “crsp” computes the R squared. The formula is 1 - esd/esd0, in this case 1 - 6.71144/7.315046 (esd is the error variance for the model and esd0 is the error variance for the empty model). Note if negative, it is replaced by zero.
The squared multiple correlation is .1582226.
To look at the measure of nonindependence, called rho, is equal to .6166651. This partial intraclass correlation is the correlations between the two members’ Satisfaction, controlling for both of their levels of Attachment Anxiety.
Need to first compute mean and standard deviation for all variables for the dataset analyzed. Then compute the standardized variables and then rerun the APIM analysis.
tempd=na.omit(riggsp[,c("Dyad","Anxiety_A","Anxiety_P","Sat_A")])
# (tempd is the dataset of cases used for the APIM run,
# all missing cases are excluded)
meanAnx = mean(tempd$Anxiety_A, na.rm=TRUE)
sdAnx = sd(tempd$Anxiety_A, na.rm=TRUE )
meanSat = mean(tempd$Sat_A, na.rm=TRUE)
sdSat = sd(tempd$Sat_A, na.rm=TRUE)
riggsp$Anxiety_As = (riggsp$Anxiety_A - meanAnx)/sdAnx
riggsp$Anxiety_Ps = (riggsp$Anxiety_P - meanAnx)/sdAnx
riggsp$Sat_As = (riggsp$Sat_A - meanSat)/sdSat
apimis = gls(Sat_As ~ Anxiety_As + Anxiety_Ps ,
na.action=na.omit,
correlation=corCompSymm (form=~1|Dyad),
data=riggsp)
smallsummary(apimis)## Correlation structure of class corCompSymm representing
## Rho
## 0.6166651
##
## Residual standard error: 0.9184858
##
## Value Std.Error t-value p-value
## (Intercept) 0.0000 0.0663 0.0000 1
## Anxiety_As -0.3044 0.0501 -6.0769 0
## Anxiety_Ps -0.2173 0.0501 -4.3376 0## 2.5 % 97.5 %
## (Intercept) -0.1300 0.1300
## Anxiety_As -0.4026 -0.2062
## Anxiety_Ps -0.3155 -0.1191Intercept = 0.00000: Is zero and should be zero because all variables are standardized..
The Actor Effect of Anxiety = -0.304: The actor beta coefficient.
The Partner Effect of Anxiety = -0.217: The partner beta coefficient.
Note that one minus the error variance (1 - .91848580x.9184858 = .1563838) closely approximates the multiple correlation squared (0.1582226)
tempd=na.omit(riggsp[,c("Dyad","Anxiety_A","Anxiety_P","Sat_A")])
# (tempd is the dataset of cases used for the APIM run,
# all missing cases are excluded)
meanAnx = mean(tempd$Anxiety_A, na.rm=TRUE)
sdAnx = sd(tempd$Anxiety_A, na.rm=TRUE )
meanSat = mean(tempd$Sat_A, na.rm=TRUE)
sdSat = sd(tempd$Sat_A, na.rm=TRUE)
tnobs=nrow(riggsp)
nobs=nrow(tempd)
tmiss=tnobs-nobs
ndyad <- length(unique(tempd$Dyad))
nsing=nobs-2*ndyadThe mean for Attachment Anxiety: 2.7796778
The mean for Satisfaction: 44.816129
The standard deviation for Attachment Anxiety: 1.2433165
The standard deviation for Satisfaction: 7.3070696
The total number of observations: 310
The total number of missing observations: 0
The total number of dyads: 155
The total number of singles: 0