library(agricolae)
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(ART)
library(ARTool)
library(car)
## Loading required package: carData
library(ggplot2)
ds=read.csv(file = "dataset.csv")
str(ds)
## 'data.frame': 60 obs. of 4 variables:
## $ Respuesta: num 0.236 0.238 0.233 0.214 0.215 ...
## $ Longitud : int 1 1 1 2 2 2 3 3 3 4 ...
## $ Masa : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Angulo : int 1 1 1 1 1 1 1 1 1 1 ...
View(ds)
f.long=factor(ds$Longitud)
f.masa=factor(ds$Masa)
f.angulo=factor(ds$Angulo)
# Gráficos de caja para determinar si existen outliers
boxplot(ds$Respuesta~(f.long+f.masa+f.angulo),col="darkgray",xlab = "Combinaciones",ylab = "Periodo (T)",cex.names=1,las=2)
g.1=boxplot(ds$Respuesta~f.long,xlab = "Factor longitud de cuerda",ylab="Periodo (T)",main="Respuesta por cada nivel")
g.2=boxplot(ds$Respuesta~f.masa,xlab = "Factor masa",ylab="Periodo (T)",main="Respuesta por cada nivel")
g.3=boxplot(ds$Respuesta~f.angulo,xlab = "Factor angulo",ylab="Periodo (T)",main="Respuesta por cada nivel")
#Análisis estadístico
modelo.pendulo=lm(ds$Respuesta~(f.long*f.masa*f.angulo))
summary(modelo.pendulo)
##
## Call:
## lm(formula = ds$Respuesta ~ (f.long * f.masa * f.angulo))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0060333 -0.0009167 -0.0000333 0.0009250 0.0069333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.2353667 0.0015614 150.740 < 2e-16 ***
## f.long2 -0.0204667 0.0022082 -9.269 1.65e-11 ***
## f.long3 -0.0394000 0.0022082 -17.843 < 2e-16 ***
## f.long4 -0.0645000 0.0022082 -29.210 < 2e-16 ***
## f.long5 -0.0907000 0.0022082 -41.075 < 2e-16 ***
## f.masa2 -0.0023667 0.0022082 -1.072 0.290
## f.angulo2 0.0008333 0.0022082 0.377 0.708
## f.long2:f.masa2 0.0017667 0.0031228 0.566 0.575
## f.long3:f.masa2 0.0021667 0.0031228 0.694 0.492
## f.long4:f.masa2 0.0009667 0.0031228 0.310 0.759
## f.long5:f.masa2 0.0018000 0.0031228 0.576 0.568
## f.long2:f.angulo2 0.0017333 0.0031228 0.555 0.582
## f.long3:f.angulo2 0.0011667 0.0031228 0.374 0.711
## f.long4:f.angulo2 0.0017000 0.0031228 0.544 0.589
## f.long5:f.angulo2 0.0007333 0.0031228 0.235 0.816
## f.masa2:f.angulo2 0.0047333 0.0031228 1.516 0.137
## f.long2:f.masa2:f.angulo2 -0.0027000 0.0044163 -0.611 0.544
## f.long3:f.masa2:f.angulo2 -0.0029667 0.0044163 -0.672 0.506
## f.long4:f.masa2:f.angulo2 -0.0037333 0.0044163 -0.845 0.403
## f.long5:f.masa2:f.angulo2 -0.0051333 0.0044163 -1.162 0.252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002704 on 40 degrees of freedom
## Multiple R-squared: 0.9953, Adjusted R-squared: 0.9931
## F-statistic: 445.7 on 19 and 40 DF, p-value: < 2.2e-16
anvar=aov(modelo.pendulo)
summary(anvar)
## Df Sum Sq Mean Sq F value Pr(>F)
## f.long 4 0.06177 0.015444 2111.504 < 2e-16 ***
## f.masa 1 0.00000 0.000000 0.026 0.871885
## f.angulo 1 0.00012 0.000119 16.232 0.000244 ***
## f.long:f.masa 4 0.00001 0.000001 0.202 0.935556
## f.long:f.angulo 4 0.00001 0.000002 0.296 0.878698
## f.masa:f.angulo 1 0.00001 0.000013 1.711 0.198351
## f.long:f.masa:f.angulo 4 0.00001 0.000003 0.362 0.833878
## Residuals 40 0.00029 0.000007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Pruebas de rango múltiple
rank.long=duncan.test(ds$Respuesta,trt=f.long,MSerror =0.000007,DFerror = 40,alpha = 0.05 )
rank.long$groups
## ds$Respuesta groups
## 1 0.2357833 a
## 2 0.2163917 b
## 3 0.1973083 c
## 4 0.1716833 d
## 5 0.1450667 e
rank.angulo=duncan.test(ds$Respuesta,trt=f.angulo,MSerror =0.000007,DFerror = 40,alpha = 0.05 )
rank.angulo$groups
## ds$Respuesta groups
## 2 0.1946533 a
## 1 0.1918400 b
plot(rank.long)
plot(rank.angulo)
shapiro.test(residuals.lm(modelo.pendulo))
##
## Shapiro-Wilk normality test
##
## data: residuals.lm(modelo.pendulo)
## W = 0.96049, p-value = 0.0497
g.norm=qqplot(seq(-3,3,by=0.01),residuals.lm(modelo.pendulo),xlab = "Residuales Estándar",ylab = "Residuales del modelo",pch="*",main="Gráfica de Probabilidad Normal")
qqline(residuals.lm(modelo.pendulo),col="red",lty=2)
abline(h=0,col="blue",lty=3)
dwtest(modelo.pendulo)
##
## Durbin-Watson test
##
## data: modelo.pendulo
## DW = 2.5099, p-value = 0.2722
## alternative hypothesis: true autocorrelation is greater than 0
residualPlot(modelo.pendulo,type="rstandard",xlab="Valores ajustados",ylab="Residuales Estándar",pch="*",main="Residuales contra ajustes")
boxplot(modelo.pendulo$residuals~f.long,xlab = "Factor longitud de cuerda",ylab="Residuales Estándar",main="Residuales contra niveles")
abline(h=0,lty=3,col="red")
boxplot(modelo.pendulo$residuals~f.masa,xlab = "Factor masa",ylab="Residuales Estándar",main="Residuales contra niveles")
abline(h=0,lty=3,col="red")
boxplot(modelo.pendulo$residuals~f.angulo,xlab = "Factor ángulo",ylab="Residuales Estándar",main="Residuales contra niveles")
abline(h=0,lty=3,col="red")
ncvTest(modelo.pendulo)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.03921503, Df = 1, p = 0.84302
#Dado que no existen condiciones de normalidad, se prueba la libreria ARtool, para realizar la transformacion de rango alineada
ds.art=data.frame(ds$Respuesta,as.factor(ds$Longitud),as.factor(ds$Masa),as.factor(ds$Angulo))
names(ds.art)=list("Respuesta","Longitud","Masa","Angulo")
View(ds.art)
modelo.art=art(Respuesta~(Longitud*Masa*Angulo),data = ds.art)
anova.art=anova(modelo.art)
print(modelo.art)
## Aligned Rank Transform of Factorial Model
##
## Call:
## art(formula = Respuesta ~ (Longitud * Masa * Angulo), data = ds.art)
##
## Column sums of aligned responses (should all be ~0):
## Longitud Masa Angulo
## 0 0 0
## Longitud:Masa Longitud:Angulo Masa:Angulo
## 0 0 0
## Longitud:Masa:Angulo
## 0
##
## F values of ANOVAs on aligned responses not of interest (should all be ~0):
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
print(anova.art)
## Analysis of Variance of Aligned Rank Transformed Data
##
## Table Type: Anova Table (Type III tests)
## Model: No Repeated Measures (lm)
## Response: art(Respuesta)
##
## Df Df.res F value Pr(>F)
## 1 Longitud 4 40 2.4615e+02 < 2.22e-16 ***
## 2 Masa 1 40 6.1064e-04 0.98041
## 3 Angulo 1 40 1.9825e+01 6.6477e-05 ***
## 4 Longitud:Masa 4 40 4.5141e-01 0.77074
## 5 Longitud:Angulo 4 40 7.2229e-01 0.58185
## 6 Masa:Angulo 1 40 2.8290e+00 0.10037
## 7 Longitud:Masa:Angulo 4 40 7.2438e-01 0.58048
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Analisis post-hoc
# Efectos principales
contr.longitud=art.con(modelo.art,"Longitud",adjust="tukey")
## NOTE: Results may be misleading due to involvement in interactions
print(contr.longitud)
## contrast estimate SE df t.ratio p.value
## Longitud1 - Longitud2 12 1.71 40 7.016 <.0001
## Longitud1 - Longitud3 24 1.71 40 14.033 <.0001
## Longitud1 - Longitud4 36 1.71 40 21.049 <.0001
## Longitud1 - Longitud5 48 1.71 40 28.066 <.0001
## Longitud2 - Longitud3 12 1.71 40 7.016 <.0001
## Longitud2 - Longitud4 24 1.71 40 14.033 <.0001
## Longitud2 - Longitud5 36 1.71 40 21.049 <.0001
## Longitud3 - Longitud4 12 1.71 40 7.016 <.0001
## Longitud3 - Longitud5 24 1.71 40 14.033 <.0001
## Longitud4 - Longitud5 12 1.71 40 7.016 <.0001
##
## Results are averaged over the levels of: Masa, Angulo
## P value adjustment: tukey method for comparing a family of 5 estimates
contr.masa=art.con(modelo.art,"Masa",adjust="tukey")
## NOTE: Results may be misleading due to involvement in interactions
print(contr.masa)
## contrast estimate SE df t.ratio p.value
## Masa1 - Masa2 0.133 5.4 40 0.025 0.9804
##
## Results are averaged over the levels of: Longitud, Angulo
contr.angulo=art.con(modelo.art,"Angulo",adjust="tukey")
## NOTE: Results may be misleading due to involvement in interactions
print(contr.angulo)
## contrast estimate SE df t.ratio p.value
## Angulo1 - Angulo2 -19.9 4.46 40 -4.453 0.0001
##
## Results are averaged over the levels of: Longitud, Masa