# R/Fract.Poly.R In CorrMixed: Estimate Correlations Between Repeatedly Measured Endpoints (E.g., Reliability) Based on Linear Mixed-Effects Models

#### Documented in Fract.Poly

```Fract.Poly <- function(Covariate, Outcome,
S=c(-2,-1,-0.5,0,0.5,1,2,3), Max.M=5, Dataset) {

if (Max.M > 5){
cat("The maximum allowed Max.M is 5. The requested Max.M is replaced by Max.M=5\n")
Max.M <- c(5)
}

Covariate <- Dataset[,paste(substitute(Covariate))]
Outcome <- Dataset[,paste(substitute(Outcome))]
n=length(Outcome)
Results_m1 <- Results_m2 <- Results_m3 <- Results_m4 <- Results_m5 <- NULL

for (j in 1: Max.M){
r <- length(S)  #aantal in S
n <- j    #aantal k
aantal <- factorial(n+r-1)/
(factorial(n) * factorial(r-1))
if (aantal > 5000){
cat("Warning: A total of ", aantal, " models of degree ", i, " can be fitted. Consider using \n", sep= "")
cat("a more restricted set S or a lower Max.M \n\n", sep= "") }
}

Covariate[Covariate==0] <- 0.00000001

#m1
if (Max.M>=1){
Results_m1 <- NULL
for (i in 1:length(S)){
if  (S[i] != 0){term1 <- Covariate ** S[i]}
if  (S[i] == 0){term1 <- log(Covariate)}
fit <- lm(Outcome~term1, data=Dataset)
AIC <- extractAIC(fit)
Results_m1 <- rbind(Results_m1, cbind(S[i], AIC))
colnames(Results_m1) <- c("power1", "AIC")
}
}

# m=2
if (Max.M>=2){
Results_m2 <- NULL
for (i in 1:length(S)){
for (j in 1:length(S)){
if  (S[i] != 0){term1 <- Covariate ** S[i]}
if  (S[i] == 0){term1 <- log(Covariate)}
if  (S[j] != 0){term2 <- Covariate ** S[j]}
if  (S[j] == 0){term2 <- log(Covariate)}
if  (S[i] == S[j]){term2 <- (Covariate ** S[j]) * log(Covariate)}
fit <-
lm(Outcome~term1+term2, data=Dataset)
AIC <- extractAIC(fit)
Results_m2 <- rbind(Results_m2, cbind(S[i], S[j], AIC))
}
colnames(Results_m2) <- c("power1", "power2", "AIC")
}
}
#m=3
if (Max.M>=3){
Results_m3 <- NULL
for (i in 1:length(S)){
for (j in 1:length(S)){
for (k in 1: length(S)){
if  (S[i] != 0){term1 <- Covariate ** S[i]}
if  (S[i] == 0){term1 <- log(Covariate)}
if  (S[j] != 0){term2 <- Covariate ** S[j]}
if  (S[j] == 0){term2 <- log(Covariate)}
if  (S[k] != 0){term3 <- Covariate ** S[k]}
if  (S[k] == 0){term3 <- log(Covariate)}
if  (S[i] == S[j]){term2 <- (Covariate ** S[j]) * log(Covariate)}
if  (S[i] == S[k]){term2 <- (Covariate ** S[k]) * log(Covariate)}
if  (S[j] == S[k]){term2 <- (Covariate ** S[k]) * log(Covariate)}
fit <- lm(Outcome~term1+term2+term3, data=Dataset)
if (NA %in% fit\$coeff == FALSE){
AIC <- extractAIC(fit)
Results_m3 <- rbind(Results_m3, cbind(S[i], S[j], S[k], AIC))}

}
}
colnames(Results_m3) <- c("power1", "power2", "power3", "AIC")
}
}

# m=4
if (Max.M>=4){
Results_m4 <- NULL
for (i in 1:length(S)){
for (j in 1:length(S)){
for (k in 1: length(S)){
for (l in 1: length(S)){
if  (S[i] != 0){term1 <- Covariate ** S[i]}
if  (S[i] == 0){term1 <- log(Covariate)}
if  (S[j] != 0){term2 <- Covariate ** S[j]}
if  (S[j] == 0){term2 <- log(Covariate)}
if  (S[k] != 0){term3 <- Covariate ** S[k]}
if  (S[k] == 0){term3 <- log(Covariate)}
if  (S[l] != 0){term4 <- Covariate ** S[l]}
if  (S[l] == 0){term4 <- log(Covariate)}
if  (S[i] == S[j]){term2 <- (Covariate ** S[j]) * log(Covariate)}
if  (S[i] == S[k]){term2 <- (Covariate ** S[k]) * log(Covariate)}
if  (S[i] == S[l]){term2 <- (Covariate ** S[l]) * log(Covariate)}
if  (S[j] == S[k]){term2 <- (Covariate ** S[k]) * log(Covariate)}
if  (S[j] == S[l]){term2 <- (Covariate ** S[l]) * log(Covariate)}
if  (S[k] == S[l]){term2 <- (Covariate ** S[l]) * log(Covariate)}
fit <- lm(Outcome~term1+term2+term3+term4, data=Dataset)

if (NA %in% fit\$coeff == FALSE){
AIC <- extractAIC(fit)
Results_m4 <- rbind(Results_m4, cbind(S[i], S[j], S[k], S[l], AIC))}
}
}
}
colnames(Results_m4) <- c("power1", "power2", "power3", "power4", "AIC")
}
}

if (Max.M>=5){
Results_m5 <- NULL
for (i in 1:length(S)){
for (j in 1:length(S)){
for (k in 1: length(S)){
for (l in 1: length(S)){
for (m in 1: length(S)){
if  (S[i] != 0){term1 <- Covariate ** S[i]}
if  (S[i] == 0){term1 <- log(Covariate)}
if  (S[j] != 0){term2 <- Covariate ** S[j]}
if  (S[j] == 0){term2 <- log(Covariate)}
if  (S[k] != 0){term3 <- Covariate ** S[k]}
if  (S[k] == 0){term3 <- log(Covariate)}
if  (S[l] != 0){term4 <- Covariate ** S[l]}
if  (S[l] == 0){term4 <- log(Covariate)}
if  (S[m] != 0){term5 <- Covariate ** S[m]}
if  (S[m] == 0){term5 <- log(Covariate)}
if  (S[i] == S[j]){term2 <- (Covariate ** S[j]) * log(Covariate)}
if  (S[i] == S[k]){term2 <- (Covariate ** S[k]) * log(Covariate)}
if  (S[i] == S[l]){term2 <- (Covariate ** S[l]) * log(Covariate)}
if  (S[i] == S[m]){term2 <- (Covariate ** S[m]) * log(Covariate)}
if  (S[j] == S[k]){term2 <- (Covariate ** S[k]) * log(Covariate)}
if  (S[j] == S[l]){term2 <- (Covariate ** S[l]) * log(Covariate)}
if  (S[j] == S[m]){term2 <- (Covariate ** S[m]) * log(Covariate)}
if  (S[k] == S[l]){term2 <- (Covariate ** S[l]) * log(Covariate)}
if  (S[k] == S[m]){term2 <- (Covariate ** S[m]) * log(Covariate)}
if  (S[l] == S[m]){term2 <- (Covariate ** S[m]) * log(Covariate)}
fit <- lm(Outcome~term1+term2+term3+term4+term5, data=Dataset)

if (NA %in% fit\$coeff == FALSE){
AIC <- extractAIC(fit)
Results_m5 <- rbind(Results_m5, cbind(S[i], S[j], S[k], S[l], S[m], AIC))}

}
}
}
}
}
colnames(Results_m5) <- c("power1", "power2", "power3", "power4", "power5", "AIC")
}

fit <- list(Results.M1=Results_m1, Results.M2=Results_m2, Results.M3=Results_m3,
Results.M4=Results_m4, Results.M5=Results_m5)

class(fit) <- "Fract.Poly"
fit
}
```

## Try the CorrMixed package in your browser

Any scripts or data that you put into this service are public.

CorrMixed documentation built on Nov. 10, 2019, 9:06 a.m.