Script Taper_editado compl

Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original

####################################################
# Ajuste de Modelos de Taper # 
####################################################
#1) Importando arquivo de dados:
T <- read.table(file="C:/User_R/dadostaper_edit.txt",head=T)
T
names(T)
#2) Ajustando modelos Taper:
di <- (T$di)
d <- (T$d)
did2 <- (di/d)^2
hih <- (T$hi/T$h)
did <- (di/d)
ldi <- log10(T$di)
ld <- log10(T$d)
lhhi <- log10(T$h-T$hi)
lh <- log10(T$h)
#a) Modelo 1 (md1): Demaerschalk (1971)
md1 <- nls(did2 ~ 10^(2*b0)*d^(2*b1-2)*(h-hi)^(2*b2)*h^(2*b3), data=T, start=list (b0=0.1, b1=0.1, b2=0.1, b3=0.1))
md1
summary(md1)
names(md1)
residuals(md1)
md1est <- (10^(0.1450)*d^(0.9114)*(T$h-T$hi)^(0.8386)*(T$h)^(-0.8356))
md1est
Emd1 <- (T$di) - md1est
Emd1
summary(Emd1)
EPmd1 <- (Emd1/(T$di))*100
EPmd1
summary(EPmd1)
mdpmd1 <- mean(EPmd1)
mdpmd1
c <- (sum(T$di))^2/282
c
sqtotal<- (sum(T$di^2))-c
sqtotal
sqresmd1<- sum(Emd1^2)
sqresmd1
glresmd1<- 278
glresmd1
n1<-282-1
r2md1<- 1-(sqresmd1/sqtotal)
r2md1
r2ajust<- 1-((1-r2md1)*(n1/glresmd1))
r2ajust
syxmd1<- sqrt(sqresmd1/glresmd1)
syxmd1
syxpmd1<- (syxmd1/mean(T$di))*100
syxpmd1
#b) Modelo 2 (md2): Kozak et. al (1969)
md2 <- lm((T$did2) ~ (T$hih)*((T$hih)^2))
md2
summary(md2)
names(md2)
residuals(md2)
md2est <- (fitted.values(md2))
md2est
Emd2 <- (T$di) - md2est
Emd2
summary(Emd2)
EPmd2 <- (Emd2/(T$di))*100
EPmd2
summary(EPmd2)
mdpmd2 <- mean(EPmd2)
mdpmd2
c <- (sum(T$di))^2/282
c
sqtotal<- (sum(T$di^2))-c
sqtotal
sqresmd2<- sum(Emd2^2)
sqresmd2
glresmd2<- 281
glresmd2
n1<-282-1
r2md2<- 1-(sqresmd2/sqtotal)
r2md2
r2ajust<- 1-((1-r2md2)*(n1/glresmd2))
r2ajust
syxmd2<- sqrt(sqresmd2/glresmd2)
syxmd2
syxpmd2<- (syxmd2/mean(T$di))*100
syxpmd2
dado<- read.table(file="C:/Users/natan/Desktop/soft R/dadostaper_edit.txt",head=T)
names (dado)
###mdp
mod2<- lm(( did2 ~ hih + hih2), data=dado)
diestm2<- sqrt((dado$di^2)*(fitted.values(mod2)))
emod2<- dado$di-diestm2
epmod2<- (emod2/dado$di)*100
mepmod2<- mean(epmod2)
mepmod2
##syx
dimedia <- mean (dado$di)
sqresm2<- sum(emod2^2)
glres2<- (280)
syxm2<- sqrt(sqresm2/glres)
syxpm2<- (syxm2/dimedia)*100
syxpm2
####r2
c <- ((sum(dado$di))^2)/282
di2<- (dado$di)^2
sqtotais<- (sum(di2)-c)
r2m2<- (1-(sqresm2/sqtotais))
# r2 ajustado
r2ajustadom2 <- (1- (1-r2m2)*(282/280))####68 é o numero de arvores - 1 , o 66 é o numero de arvores menos o numero d parametros.
########################################gráficos 
par(mfrow=c(1,2))
plot (diestm2,emod2,xlab="di (estimado)",ylab="erro absoluto(m)",pch=18,main="modelo kozak et al",xlim=c(0,30),ylim=c(-10,35))
plot (dado$d,epmod2,xlab="dap",ylab="erro percentual",pch=18,main="modelo kozak et al",xlim=c(0,30),ylim=c(-20,70))
###############################m3 
names (dado)
did<- dado$d/dado$di
###mdp
mod3<- lm(( did ~ hih + hih2 + hih3+hih4+hih5), data=dado )
diestm3<- (dado$di*(fitted.values(mod3)))
emod3<- (dado$di - diestm3)
epmod3<- (emod3/dado$di)*100
mepmod3<- mean(epmod3)
mepmod3
##syx
glres3<- (282-5)
sqresm3<- sum(emod3^2)
syxm3<- sqrt(sqresm3/glres3)
syxpm3<- (syxm3/dimedia)*100
syxpm3
####r2
r2m3<- (1-(sqresm3/sqtotais))
# r2 ajustado
r2ajustadom3 <- (1- (1-r2m3)*(282/277))####68 é o numero de arvores - 1 , o 66 é o numero de arvores menos o numero d parametros.
r2ajustadom3
par(mfrow=c(1,2))
plot (diestm3,emod3,xlab="di (estimado)",ylab="erro absoluto(m)",pch=18,main="polinomio do 5 grau",xlim=c(-10,30),ylim=c(-30,35))
plot (dado$d,epmod3,xlab="dap",ylab="erro percentual",pch=18,main="polinomio do 5 grau",xlim=c(0,30),ylim=c(-500,35))
#############################################m4
names (dado)
###mdp
mod4<- lm(( logdi ~ logd + loghhi +logh), data=dado )
diestm4<- (exp(fitted.values(mod4)))
emod4<- (dado$di - diestm4)
epmod4<- (emod4/dado$di)*100
mepmod4<- mean(epmod4)
mepmod4
##syx
glres4<- (279)
sqresm4<- sum(emod4^2)
syxm4<- sqrt(sqresm4/glres)
syxpm4<- (syxm4/dimedia)*100
syxpm4
####r2
r2m4<- (1-(sqresm4/sqtotais))
# r2 ajustado
r2ajustadom4 <- (1- (1-r2m4)*(282/279))####68 é o numero de arvores - 1 , o 66 é o numero de arvores menos o numero d parametros.
r2ajustadom4
par(mfrow=c(1,2))
plot (diestm4,emod4,xlab="di (estimado)",ylab="erro absoluto(m)",pch=18,main="omerad",xlim=c(0,30),ylim=c(0,35))
plot (dado$d,epmod4,xlab="dap",ylab="erro percentual",pch=18,main="omerad",xlim=c(0,50),ylim=c(40,100))
###############################m5
names (dado)
###mdp
mod5<- nls(( did2 ~ (((h-hi)/(h-1.3))^2*b1)), data=dado , start=list(b0=0.01, b1=0.1, b2=0.1, b3=0.1))
diestm5<- sqrt((dado$di^2)*(fitted.values(mod5))
emod5<- (dado$di - diestm4)
epmod4<- (emod4/dado$di)*100
mepmod4<- mean(epmod4)
mepmod3
##syx
glres4<- (279)
sqresm4<- sum(emod4^2)
syxm4<- sqrt(sqresm4/glres)
syxpm4<- (syxm4/dimedia)*100
syxpm4
####r2
r2m4<- (1-(sqresm4/sqtotais))
# r2 ajustado
r2ajustadom4 <- (1- (1-r2m4)*(282/279))####68 é o numero de arvores - 1 , o 66 é o numero de arvores menos o numero d parametros.
r2ajustadom4