AJUSTES DE MODELOS POR CURVA GUIA
Curso recivido por el Dr Juan Gabriel. IPS-USC-LUGO, ESPAÑA, OCTUBRE 2015
MANUAL DE AJUSTE DE MODELOS (PASO POR PASO)
caragar datos en eviroment/importar dataset/poner un nombre y yes en cabecera/importar
Datos <- read.csv("C:/LUGO/CURSOR/Bases_datos/TreeVolumes.csv")
View(Datos)#ver los datos
View(Datos$volume.dm3)#Ver la columna volumen unicamnete
#agregar columas
Datos$d2<-Datos$dbh.cm**2 #Calcular la columna d2
Datos$vd2<-Datos$volume.dm3*Datos$d2 #Calcular la columna vd2
Datos$d4<-Datos$dbh**4 #Calcular la columna d4
# Calcualr el valor b
mumb<-sum(Datos$vd2)-sum(Datos$volume.dm3)*sum(Datos$d2)/length(Datos$d2)
denb<-sum(Datos$d4)-(sum(Datos$d2))**2/length(Datos$d2)
b<-mumb/denb
#Calcular el valor a
a<-mean(Datos$volume.dm3)-b*mean(Datos$d2)
#por lo tanto tendremos la recta V=-115.1388+0.6862*d2 en cm.AJUSTE EN R DE UN MODELO LINEAL
variables dependientes(Y)=variables independientes(x+z)
recta<-lm(volume.dm3~d2, Datos) #Es igual lm(Datos$volume.dm3~Datos$d2), se recomienda la primera obción.Al poner nombre(recta) aparece en Evironmente los datos con una serie de descripción. por ejemplo en el coeficientes nos muetra el valor de los coeficientes a y b. siguiente scrib.
recta$coefficients #para que nos muestre los coeficientes de la recta.## (Intercept) d2
## -115.1388740 0.6862423summary(recta) #para que nos muestre el resumen del modelo.##
## Call:
## lm(formula = volume.dm3 ~ d2, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -99.011 -51.517 -0.093 51.615 131.836
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -115.13887 30.69128 -3.752 0.000815 ***
## d2 0.68624 0.04549 15.086 5.66e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 65.74 on 28 degrees of freedom
## Multiple R-squared: 0.8904, Adjusted R-squared: 0.8865
## F-statistic: 227.6 on 1 and 28 DF, p-value: 5.658e-15en esta parte nos muestra el residual mas pequeño, el mas grande, primer cuartil, mediana, 3q. El residual standartd error es llamada tambien root mean squarer error o la raiz del error medio cuadratico(RMSE). El multiple R-squared nos muestra el error, por lo tanto el error seria lo faltante(10.06)#?revisar el adjusted r-squared es la que se utiliza para comparar modelos, ya que toma en cuenta los variables empleados.
OBTENER LOS VARIABLES MANUALES
obtener el r2, los residuales se pueden ver en la caja (lado derecho).
numr2<-sum((recta$residuals)**2)
denr2<-sum((Datos$volume.dm3-mean(Datos$volume.dm3))**2)
r2<-1-numr2/denr2
#r ajustado
numr2a<-numr2/28
denr2a<-denr2/29
r2a<-1-numr2a/denr2a
rmse<-sqrt(numr2/28)
anova(recta)## Analysis of Variance Table
##
## Response: volume.dm3
## Df Sum Sq Mean Sq F value Pr(>F)
## d2 1 983438 983438 227.58 5.658e-15 ***
## Residuals 28 120994 4321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sqrt(4321)## [1] 65.73431AJUSTAR MODELO v=a+bd2+ctree
recta2<-lm(volume.dm3~d2+tree, Datos)
summary(recta2)##
## Call:
## lm(formula = volume.dm3 ~ d2 + tree, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -96.211 -63.983 -3.883 39.874 126.608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -132.83369 39.15559 -3.392 0.00215 **
## d2 0.68892 0.04601 14.974 1.34e-14 ***
## tree 1.03433 1.40241 0.738 0.46716
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 66.28 on 27 degrees of freedom
## Multiple R-squared: 0.8926, Adjusted R-squared: 0.8847
## F-statistic: 112.2 on 2 and 27 DF, p-value: 8.279e-14coparar los l?os modelos
modelo 1 r2=0.8904, r2 adj=0.8865, Residual standard error (RMSE): 65.74
modelo 2 r2=0.8926,, r2 adj=0.8847,
este modelo (2)es mejor segun r, pero r2 ajustado nos dice que el primero (1) es mejor, por lo tanto el mejor es el primero (1).
quitar el termino independiente y compara modelos.
recta3<-lm(volume.dm3~d2-1,Datos)
summary(recta3)##
## Call:
## lm(formula = volume.dm3 ~ d2 - 1, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -126.53 -74.25 -40.17 42.55 147.88
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## d2 0.52918 0.02143 24.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 79.18 on 29 degrees of freedom
## Multiple R-squared: 0.9546, Adjusted R-squared: 0.9531
## F-statistic: 610 on 1 and 29 DF, p-value: < 2.2e-16modelo 3 r2=0.9546,, r2 adj=0.9531, Residual standard error (RMSE): 76.18
Este modelo nos muestra ser la mejor sun embargo el RMSE mustra que el error es mayor por lo tanto el mejor es la primera por el RMSE menor.
numr22<-sum((recta3$residuals)**2)
denr22<-sum((Datos$volume.dm3-mean(Datos$volume.dm3))**2)
r22<-1-numr22/denr22AJSTAR MODELO v=a+bd2h
Datos$d2h<-Datos$d2*Datos$height.m
recta4<-lm(volume.dm3~d2h, Datos)
summary(recta4)##
## Call:
## lm(formula = volume.dm3 ~ d2h, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -71.689 -10.537 -0.636 9.910 56.056
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.426e+01 9.663e+00 1.476 0.151
## d2h 3.596e-02 9.862e-04 36.459 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.53 on 28 degrees of freedom
## Multiple R-squared: 0.9794, Adjusted R-squared: 0.9786
## F-statistic: 1329 on 1 and 28 DF, p-value: < 2.2e-16modelo 4 r2=0.9794, r2 adj=0.9794, Residual standard error (RMSE): 28.53
este modelo es el mejor
Cosas que hay que mirar siempre en modelos lineales
estra trabajando con los mismos datos
que todos tengan un termino independiente
que la variable dependiente se ala misma(v)
quitar el temino independiete, por eso esta el -1
recta5<-lm(volume.dm3~d2h-1, Datos)
summary(recta5)##
## Call:
## lm(formula = volume.dm3 ~ d2h - 1, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -73.851 -5.979 6.771 18.455 52.602
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## d2h 0.0371823 0.0005422 68.57 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 29.1 on 29 degrees of freedom
## Multiple R-squared: 0.9939, Adjusted R-squared: 0.9937
## F-statistic: 4702 on 1 and 29 DF, p-value: < 2.2e-16esta pero por no tener variable independiente
calcular el r2 verdadero del modelo
OBTENER LOS VARIABLES MANUALES
obtener el r2, los residuales enla caja (lado derecho) muestran los errores
numr5<-sum((recta5$residuals)**2)
denr5<-sum((Datos$volume.dm3-mean(Datos$volume.dm3))**2)
r2<-1-numr5/denr5r ajustado
numr5a<-numr5/28 #(20-(n-p)),n=numero total de observaciones-p(2)
denr5a<-denr5/29 #
r5a<-1-numr5a/denr5a
rmse<-sqrt(numr5/28)
anova(recta5)## Analysis of Variance Table
##
## Response: volume.dm3
## Df Sum Sq Mean Sq F value Pr(>F)
## d2h 1 3981455 3981455 4702 < 2.2e-16 ***
## Residuals 29 24556 847
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sqrt(877) #Comprobar el rmse que coincidan(sqrt(el valor de numr5a))## [1] 29.61419DIA 1/10/15
ajustar un modelo no lineal v=ad(b)h(c) con nls
Datos <- read.csv("C:/LUGO/CURSOR/Bases_datos/TreeVolumes.csv")
ajuste.nl<-nls(volume.dm3~a*dbh.cm**b*height.m**c, Datos,start=c(a=0.04, b=2, c=1))start=es el valor de partida de los parametros, dede el valor de b y c se toman de la formula general para calcular el volumen de un arbol v=s(pi/4)d2h
summary(ajuste.nl)##
## Formula: volume.dm3 ~ a * dbh.cm^b * height.m^c
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.07684 0.03096 2.482 0.0196 *
## b 1.73748 0.14275 12.171 1.79e-12 ***
## c 1.05803 0.09499 11.138 1.34e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.34 on 27 degrees of freedom
##
## Number of iterations to convergence: 7
## Achieved convergence tolerance: 4.581e-07calcular r2
numr<-sum((resid(ajuste.nl)**2)) #
sqrt(numr/27)#calculo de numr(anterior)## [1] 28.33802denr<-sum((Datos$volume.dm3-mean(Datos$volume.dm3))**2)
r2<-1-numr/denrmodelo r2=0.9803, Residual standard error (RMSE): 28.34
r2 llamarlo eficiencia de modelo (ME)
hacer un ajuste lineal
r2lineal<-lm(Datos$volume.dm3~fitted(ajuste.nl))
summary(r2lineal)#este ajuste nos da mejro ajuste que el anterior, (utilizar cualquiera de las dos)##
## Call:
## lm(formula = Datos$volume.dm3 ~ fitted(ajuste.nl))
##
## Residuals:
## Min 1Q Median 3Q Max
## -65.137 -11.306 -0.071 10.336 56.394
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.60370 9.63948 0.374 0.711
## fitted(ajuste.nl) 0.99159 0.02645 37.488 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 27.76 on 28 degrees of freedom
## Multiple R-squared: 0.9805, Adjusted R-squared: 0.9798
## F-statistic: 1405 on 1 and 28 DF, p-value: < 2.2e-16nos muestra que el a es cero(.71, significativo) y b1(.99)
si no se sabe los valores de los parametros se realiza la linealizacion del modelo.
primero que hay que comvertirlo en un modelo no lineal e intentar linealizarlo para tener solo una solucion, para el caso del modelo trabajo aplicar log.
v=ad(b)h(c)==logv=loga+blogd+clogh, fijar que cuando se tiene en la base datos un valor negativo el log no sirve poque no hay log en numeros negativos.
ocupar el mismo logaritmo para todo el modelo
Datos$logy<-log(Datos$volume.dm3)
Datos$logd<-log(Datos$dbh.cm)
Datos$logh<-log(Datos$height.m)
ajuste.log<-lm(logy~logd+logh,Datos)
summary(ajuste.log)##
## Call:
## lm(formula = logy ~ logd + logh, data = Datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.154020 -0.072066 -0.001231 0.055804 0.134780
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.66057 0.23557 -7.049 1.41e-07 ***
## logd 1.37648 0.12994 10.593 4.06e-11 ***
## logh 1.16540 0.09097 12.810 5.49e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08461 on 27 degrees of freedom
## Multiple R-squared: 0.99, Adjusted R-squared: 0.9892
## F-statistic: 1334 on 2 and 27 DF, p-value: < 2.2e-16logv=loga+blogd+clogh, como no cambiaron los coeficinetes d y h , los valores no cambian pero el coeficinete a si cambia por lo tanto hay que convertirlo (log=e^el)
exp(-1.66057)## [1] 0.1900306Estos resultados de rs no se pueden comparar con los ajustes anteriore, ya que los variables no son iguales, aqui se trabaja con log y en las anteriores no se aplica log
ajuste.nl<-nls(volume.dm3~a*dbh.cm**b*height.m**c, Datos,start=c(a=1.37648, b=1.16540, c=0.1900306))
summary(ajuste.nl)##
## Formula: volume.dm3 ~ a * dbh.cm^b * height.m^c
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.07684 0.03096 2.482 0.0196 *
## b 1.73748 0.14275 12.171 1.79e-12 ***
## c 1.05803 0.09499 11.138 1.34e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 28.34 on 27 degrees of freedom
##
## Number of iterations to convergence: 10
## Achieved convergence tolerance: 3.12e-07esto nos da lo mismo si no se linealiza
CURVAS DE CALIDAD DE ESTACION
Ley Eichorn, vio que la existencia de produccion depende de la altura media del rodal+edad o diametro cuadratico(dg), sin embargo la altura es afecta por los tratamientos por lo tanto se emplea la altura dominantes, considerando que los tratamientos no afecten los arboles dominantes, entonces produccion=altura dominante+edad.
METODOLOGIA PARA CONSTRUIR UNA CURVA
metodo de la curva guia
metodo de las curvas de diferencias algebraicas (ADA)
- metodo de las curvas de diferencias algebraicas generalizadas (GADA)-mas usada
metodo de estimacion de parametros(muy bueno pero poco utilizado)
para generar las curvas hay que considerar:
cubrir todas las edades
todos lo sitios desde los mejores hasta los peores
tener datos representativos
se pueden trabajar con:
parcelas temporales (una medicion)
parcelas pemanentes (3 o mas mediciones)
parcelas de intervalos (se miden 2 veces, coindicen con los tratamientos)
análisis de troco (el Árbol tiene que ser dominante, sin embargo no sabes de su pasado, pudo no ser el denominate y empezo a serlo despues de la corta de los dominates)
CARACTERISTICAS DE LAS CURVAS
deben de ser crecientes
tener dos puntos de inflexion (al inicio y al final).
la curvas deben de tener asintota (un punto maximo de crecimiento).
pasar por el punto 0,0 (comportamiento biologico razonable).
deben de tener un edad de referencia que se le conoce como indice de sitio (turno)
METODO DE LA CURVA GUIA
En este metodo se puede ocupar todo tipo de datos, da igual de donde venga sera dificil encontrar masas abanzadas con edades mas altas, ya que se cortan antes.
se debe elieguir un modelo que cumple con las 4 caracteristicas y por regresion no lineal,se agusta la curva guia a todos esos puntos
figar una ead de referencia, se debe escoger de la mitad del turno hacia adelante porque al inicio las plantas pueden tener lento o rapido crecimeinto po el calidad de sitio y despues se muestran el verdadero comportamiento de las cuvas guias.
ver los valores maximos y minimos a la edad de referencia, un poco por encima y por debajo del ultimo punto tenindo dos curvas(una arriba y otra abajo) en la diferencia entre ambos puntos se sacan el resto de los puntos.
caracteristicas de las curbas mas caracteristicas de la familia
tienen que ser polimorficas
tener difeentes asintontas (forma de s) (asintoto= edad maxima que alcanzan los arboles).
curvas anamorficas (proporcionales)
curvas polimorficas (no proporcionales)
las curvas son proporcionales cuando todos los factores son iguales excepto un factor, lo logico es que no sean proporcionaes
CONSTRUCCION DE CURVAS GUIA
datoscg <- read.csv("C:/LUGO/CURSOR/Bases_datos/StemAnalysis161PinusRadiata.csv")
View(datoscg)
tail(datoscg)#para que nos mueste los ultimos datos, cola.## plot tree age.yr height.m
## 1905 BAC2 161 8 6.69
## 1906 BAC2 161 9 8.17
## 1907 BAC2 161 10 9.40
## 1908 BAC2 161 11 10.40
## 1909 BAC2 161 12 11.44
## 1910 BAC2 161 13 12.89datoscg[1910, ] #para que nos muestre la ultima fila de datos## plot tree age.yr height.m
## 1910 BAC2 161 13 12.89plot(datoscg$age.yr,datoscg$height.m)#primero crear un plot para ver como se distribuyen los puntos
notar que a edades avanzadas no hay muchos datos de pinus radiata
trabajar con el modelo de HOSSFELD, uno de los de amediado del siglo XIX,
es facil de trabjar
HO=t2/(a+bt)2
datoscg$y<-datoscg$age.yr/(datoscg$height.m)**0.5 #crear una columna de Y, 0.5 es para la raiz
plot(datoscg$age.yr,datoscg$y) #crear un plot con el nuevo columa creado
#LINEALIZAR MODELO (lm)
HO=t2/(a+bt)2
raiz(HO)=t/(a+bt)
despejar a+bt=y
a+bt=t/raiz(Ho)
y=t/raiz(H0)#modelo linealizado a ajustar
reg.lineal<-lm(y~age.yr, datoscg ) #ajustar modelo
summary(reg.lineal) #el sumary es significativo y el r2 no nos sirve ya que el modelo no es lineal##
## Call:
## lm(formula = y ~ age.yr, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.14581 -0.33954 -0.07457 0.20745 2.83346
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.859381 0.021693 85.71 <2e-16 ***
## age.yr 0.140788 0.001281 109.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4803 on 1908 degrees of freedom
## Multiple R-squared: 0.8635, Adjusted R-squared: 0.8634
## F-statistic: 1.207e+04 on 1 and 1908 DF, p-value: < 2.2e-16obtener el asintota:
HO=t2/(a+bt)2
lim(t>-inf) t2/(a+bt)2
inf/inf=interminada
t2/(a+b2t^2abt)
1/(a2/t2+b^2+2ab/t)
1/(0+b^2+0) este se empleara en la siguinete parte para obtener el asintota
1/0.1408**2 #asisntota qe se va a utilizar en ln.. 0.1408 se obtiene del modelo anterior## [1] 50.44228con este ajuste nos da a=1.8593 y b=0.1408, las cuales se van a utilizar para el ajuste del nls
AJUSTAR MODELO NO LINEAL (nls)
HO=t2/(a+bt)2
reg.nl<-nls(height.m~(age.yr**2)/(a+b*age.yr)**2,datoscg, start=c(a=1.86, b=0.14))
summary(reg.nl)##
## Formula: height.m ~ (age.yr^2)/(a + b * age.yr)^2
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 1.778673 0.023522 75.62 <2e-16 ***
## b 0.141209 0.001182 119.51 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.777 on 1908 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.923e-061/0.1412**2 #asintota de b en nls## [1] 50.15689HO=t2/(a+bt)2
HO=t2(años)/(1.7787+0.1412t(años))2, con un RMSE=2.777 m
crear el resto de curva guia
fijar una edad de referencia, en este caso sera 20 años
buscar losdos puntos guia de los extremos a 20 años, seran 8m y 28 m
crear otros tres dentro del intervalo de as dos, seran a 13m, 18m y 23m
FAMILIA 1 (POLIMORFICA Y CON MULTIPLES ASINTOTAS)
HO=t2/(a+bt)2
8=202/(a+b(20))2 para la primera curva, mismo procedimiento para el resto.
utilizar el valor de a obtenido en el ajuste de nls
8=202/(1.7787+b8(20))2 sustituir en a
20/raiz(8)=1.7787+b8(20)
b8=((20/raiz(8))-1.7787)/20 #la cual se aplicara a cada recta
obtenr los dos grupos de familia
1-. a=1.778(fijo) y b=cambia, Familia polimorfica con multiples asintotas, a utilizar
2-. a=cambia y b=0.1412(fijo), polimorficas con asintota unica
b8<-(20/sqrt(8)-1.7787)/20
b13<-(20/sqrt(13)-1.7787)/20
b18<-(20/sqrt(18)-1.7787)/20
b23<-(20/sqrt(23)-1.7787)/20
b28<-(20/sqrt(28)-1.7787)/20DIBUJAR CURVAS SOBRE EL PLOT FAMILIA
x<-c(0,1:40) #hasta la edad maxima en xHO=t2/(a+bt)2 sustituir en la formula original
y8<-x**2/(1.7787+0.2646*x)**2 #en el paraemtro a= valor obtenido al ajustar,b=obtenido derivar b ventana a lado derecho
y13<-x**2/(1.7787+0.1884*x)**2
y18<-x**2/(1.7787+0.1467*x)**2
y23<-x**2/(1.7787+0.1195*x)**2
y28<-x**2/(1.7787+0.1000*x)**2
plot(datoscg$age.yr,datoscg$height.m)
lines(x,y8, col="red", lwd=2) # curvas sobre plot de distribucion
lines(x,y13, col="red",lwd=2)
lines(x,y18, col="red",lwd=2)
lines(x,y23, col="red",lwd=2)
lines(x,y28, col="red",lwd=2) 
esta familia de curva no se agusta muy bien a los puntos, en la parte de los 10 primeros años y en la quinta parte.
FAMILIA 2 (FAMILIA CON ASINTOTA UNICA)
8=202/(a+b(20))2 para la primera curva, mismo procedimiento para el resto.
8=202/(a+0.1412(20))2
8=202/(a+2.824)2
20/raiz(8)=a+2.824
a8=20/raiz(8)-2.824
a8<-(20/sqrt(8)-2.824)
a13<-(20/sqrt(13)-2.824)
a18<-(20/sqrt(18)-2.824)
a23<-(20/sqrt(23)-2.824)
a28<-(20/sqrt(28)-2.824)DIBUJAR CURVAS SOBRE EL PLOT FAMILIA
HO=t2/(a+bt)2 sustituir en la formula original
z8<-x**2/(4.2470+0.1412*x)**2 #en el paraemtro b= valor obtenido al ajustar, b=obtenido derivar a
z13<-x**2/(2.7230+0.1412*x)**2
z18<-x**2/(1.8900+0.1412*x)**2
z23<-x**2/(1.3462+0.1412*x)**2
z28<-x**2/(0.9556+0.1412*x)**2
plot(datoscg$age.yr,datoscg$height.m)
lines(x,z8, col="blue",lwd=2) ## curvas sobre plot de distribucion
lines(x,z13, col="blue",lwd=2)
lines(x,z18, col="blue", lwd=2)
lines(x,z23, col="blue",lwd=2)
lines(x,z28, col="blue",lwd=2)
esta curva se ajusta mejor a los datos.
en ambas familias se ajusto unicamente la curva central y de ahi se distrbuyen los otros.
AJUSTE DEL MODELO CHADMAN-RICHARDS NO LINEA, NO SE PUEDE LINEALIZAR
Ho=a(1-e-bt)c
Primero buscar donde esta la asintota si esta se desconoce o si se conoce de otros trabajos utilizarlos
asintota= valor al que tiende la curva cuando la edad es muy grande
asintota lim(t>inf) a(1-e-bt)c, donde:
-bt=inf e^-bt=0 1-e-bt=1c por lo tanto tendriamos
a*1, entonces nuestra asintota tiene un valor a(valor del parametro) cuando b>0
tenemos que quitar el exponente en nuestra formula original obteniendo
log Ho= log a + c log(1-e^(-bt)) donde:
log Ho=y log a=a log(1-e^(-bt))=x,
obteniedo y=a+cx, sin embargo, este modelo no es lineal debido a que no se conoce y no se puede calcular el valor ede b ya que es un parametro
datoscg$y<-log(datoscg$height.m)entonces buscaremos el mayor r2 y el menor RMSE sustituyendo valores en el parametro b hasta encontra el mejor, como se muestra acontinuación.
datoscg$x8<-log(1-exp(-8*datoscg$age.yr))#cuando el valor es 8
ajuste<-lm(y~x8, datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x8, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9556 -0.4343 0.1526 0.5909 1.1878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.362e+00 1.704e-02 138.632 < 2e-16 ***
## x8 7.092e+03 1.110e+03 6.391 2.06e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7438 on 1908 degrees of freedom
## Multiple R-squared: 0.02096, Adjusted R-squared: 0.02045
## F-statistic: 40.85 on 1 and 1908 DF, p-value: 2.062e-10datoscg$x30<-log(1-exp(-30*datoscg$age.yr))#cuando el valor es 30, donde el r2 es menor, entonces se descarta los que son >30 ya que buscamos el valor mas grande
ajuste<-lm(y~x30, datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x30, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9564 -0.4343 0.1527 0.5909 1.1878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.362e+00 1.704e-02 138.627 < 2e-16 ***
## x30 2.538e+13 3.978e+12 6.381 2.2e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7438 on 1908 degrees of freedom
## Multiple R-squared: 0.0209, Adjusted R-squared: 0.02038
## F-statistic: 40.72 on 1 and 1908 DF, p-value: 2.2e-10datoscg$x6<-log(1-exp(-6*datoscg$age.yr))# obtenermos valor > al 8 por lo tanto descartamos el valor de 8
ajuste<-lm(y~x6, datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x6, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9505 -0.4344 0.1525 0.5907 1.1877
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.36195 0.01703 138.665 < 2e-16 ***
## x6 968.30040 149.97208 6.457 1.36e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7436 on 1908 degrees of freedom
## Multiple R-squared: 0.02138, Adjusted R-squared: 0.02087
## F-statistic: 41.69 on 1 and 1908 DF, p-value: 1.355e-10datoscg$x005<-log(1-exp(-0.05*datoscg$age.yr))# obtenemos el mejor r2= .8837y rmse=0.2563(revisar)
ajuste<-lm(y~x005, datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x005, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33589 -0.12545 0.03531 0.16724 1.09789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.57249 0.01168 306.0 <2e-16 ***
## x005 1.44194 0.01197 120.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2563 on 1908 degrees of freedom
## Multiple R-squared: 0.8837, Adjusted R-squared: 0.8837
## F-statistic: 1.45e+04 on 1 and 1908 DF, p-value: < 2.2e-16con esto tenemos b=0.05, a=1.5, a=e^3.6
exp(3.6)#calcular la asintota, la cual nos da 36 y es muy alta a la recta que estamos ajustando, recordad que se está ajustando la curva central, ver en plot## [1] 36.59823plot(datoscg$age.yr, datoscg$height.m)
ajuste.nl<-nls(height.m~a*(1-exp(-b*age.yr))**c,datoscg, start=c(a=36, b=0.05, c=1.5))# ajustar el modelo no lineal, en la cual nos da una asintota de 31( valor de a)
summary(ajuste.nl)##
## Formula: height.m ~ a * (1 - exp(-b * age.yr))^c
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 31.012002 0.712315 43.54 <2e-16 ***
## b 0.069592 0.003904 17.83 <2e-16 ***
## c 1.686666 0.067585 24.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.774 on 1907 degrees of freedom
##
## Number of iterations to convergence: 4
## Achieved convergence tolerance: 9.714e-06con esto obtenemos Ho=31.012*(1-e-0.070t)1.687
por lo tanto tenemos tres obciones de familia (a,b y c)
formulas de las obciones
x<-c(0:40) #Crear un vector de acurdo a la edad mayorFAMILIA A (OBCION A)
a=variable, b=fija c=fija —–multiples asintotas
HoI/HoII=aI(1-e(-bt))c/aII(1-e-bt)c #cambios en aI y aII
HoI/HoII=aI/aII—–Curvas anamorficas o proporcionales
Del modelo origina despejamos a
a=Ho/(1-e-bt)c
a8=8/(1-e-0.070(20))1.687
a8<-8/(1-exp(-0.070*20))**1.687
a13<-13/(1-exp(-0.070*20))**1.687
a18<-18/(1-exp(-0.070*20))**1.687
a23<-23/(1-exp(-0.070*20))**1.687
a28<-28/(1-exp(-0.070*20))**1.687
y8<-a8*(1-exp(-0.070*x))**1.687
y13<-a13*(1-exp(-0.070*x))**1.687
y18<-a18*(1-exp(-0.070*x))**1.687
y23<-a23*(1-exp(-0.070*x))**1.687
y28<-a28*(1-exp(-0.070*x))**1.687
#crear curvas en el plot edad y altura
plot(datoscg$age.yr, datoscg$height.m)
lines(x, y8, col="red")
lines(x, y13, col="red")
lines(x, y18, col="red")
lines(x, y23, col="red")
lines(x, y28, col="red")
estas curvas en la parte de abajo no van muy bien
FAMILIA B (obcion b)
a=fijo, b=variable, c=fijo —–asintotas fija
HoI/HoII=a(1-e(-bIt))c/a(1-e-bIIt)c #cambios en bI y bII
HoI/HoII=[(1-e(-bIt))/(1-e(-bIIt))]^c curvas polimorficas
despejamos b del modelo original
a=Ho/(1-e-bt)c
Ho/a=(1-e-bt)c
(Ho/a)(1/c)=(1-e-bt)
(e(-bt))-t=(Ho/a)(1/c) donde e^(-bt)=-btloge y loge=1 por lo tanto
-bt=log(1-(Ho/a)^(1/c)) donde el valor de (Ho/a)^(1/c) tiene que ser menor a 1, de lo contrario nos daria error por el log.
-b=log(1-(Ho/a)^(1/2))/t
b=-log(1-(Ho/a)^(1/2))/t
crear curvas
b8<--log(1-(8/31.012)**(1/1.687))/20
b13<--log(1-(13/31.012)**(1/1.687))/20
b18<--log(1-(18/31.012)**(1/1.687))/20 # este valor nos tiene que dar un valor al rededor de 0.07
b23<--log(1-(23/31.012)**(1/1.687))/20
b28<--log(1-(28/31.012)**(1/1.687))/20
y8<-31.012*(1-exp(-b8*x))**1.687
y13<-31.012*(1-exp(-b13*x))**1.687
y18<-31.012*(1-exp(-b18*x))**1.687
y23<-31.012*(1-exp(-b23*x))**1.687
y28<-31.012*(1-exp(-b28*x))**1.687
plot(datoscg$age.yr, datoscg$height.m)
lines(x, y8, col="green")
lines(x, y13, col="green")
lines(x, y18, col="green")
lines(x, y23, col="green")
lines(x, y28, col="green")
# FAMILIA C
a=fijo, b=fijo, c=variable —–asintotas fija
HoI/HoII=a(1-e(-bt))cI/a(1-e-bIt)c1 #cambios en cI y cII
HoI/HoII=(1-e(-bt))(cI-cII) curvas polimorficas
despejamos c del modelo original
a=Ho/(1-e-bt)c despejar C
Ho/a=(1-e-bt)c
log(Ho/a)=(1-e^-bt)
c=(log (Ho/a)/log(1-e^-bt)) utilizar para sustituir
CREAR CURVAS
c8<-log(8/31.012)/log(1-exp(-0.070*20))
c13<-log(13/31.012)/log(1-exp(-0.070*20))
c18<-log(18/31.012)/log(1-exp(-0.070*20))
c23<-log(23/31.012)/log(1-exp(-0.070*20))
c28<-log(28/31.012)/log(1-exp(-0.070*20))
y8<-31.012*(1-exp(-0.069*x))**c8
y13<-31.012*(1-exp(-0.069*x))**c13
y18<-31.012*(1-exp(-0.069*x))**c18
y23<-31.012*(1-exp(-0.069*x))**c23
y28<-31.012*(1-exp(-0.069*x))**c28
plot(datoscg$age.yr, datoscg$height.m)
lines(x, y8, col="green")
lines(x, y13, col="green")
lines(x, y18, col="green")
lines(x, y23, col="green")
lines(x, y28, col="green")
DE LOS TRES FAMILIAS EL QUE MEJOR SE ADAPTA A NUESTROS PUNTOS ES LA SEGUNDA FAMILIA.
NOTA: “ES IMPOSIBLE OBTENER CURVAS POLIMORFICAS CON MULTIPLES ASINTOTAS EMPLEANDO EL METODO DE LA CURVA GUIA Y QUE TENGAN SENTIDO BIOLOGICO”, son muy pocos en las que se puede.
MODELO DE KORFT Ho=e(a+(b/(tc))
tambien conocido como modelo de Korft, bailey y clutter y/o schumacher, es uno de los modelos mas usados y sobre todo en plantas de crecimiento muy rapido
datoscg <- read.csv("C:/LUGO/CURSOR/Bases_datos/StemAnalysis161PinusRadiata.csv")Ho=e(a+(b/(tc))
log Ho= a+(a/(t^c))
Y= a+b(1/(t^c))
datoscg$y<-log(datoscg$height.m)buscad los valores mas aceptables para c
datoscg$x08<-1/datoscg$age.yr**0.8
ajuste<-lm(y~x08,datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x08, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5241 -0.1792 0.0436 0.2136 3.3866
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.38378 0.01483 228.13 <2e-16 ***
## x08 -6.45561 0.07831 -82.44 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3519 on 1908 degrees of freedom
## Multiple R-squared: 0.7808, Adjusted R-squared: 0.7807
## F-statistic: 6797 on 1 and 1908 DF, p-value: < 2.2e-16datoscg$x2<-1/datoscg$age.yr**2
ajuste<-lm(y~x2,datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x2, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7494 -0.4245 0.0682 0.4768 4.9811
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.49791 0.01528 163.45 <2e-16 ***
## x2 -7.16418 0.25364 -28.25 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6312 on 1908 degrees of freedom
## Multiple R-squared: 0.2949, Adjusted R-squared: 0.2945
## F-statistic: 797.8 on 1 and 1908 DF, p-value: < 2.2e-16datoscg$x01<-1/datoscg$age.yr**0.1 # es la mas aceptable por los valores en r2 y RMSE
ajuste<-lm(y~x01,datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x01, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.30695 -0.13530 0.03017 0.17143 0.89014
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.84084 0.08919 144.0 <2e-16 ***
## x01 -13.41617 0.11388 -117.8 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2613 on 1908 degrees of freedom
## Multiple R-squared: 0.8791, Adjusted R-squared: 0.8791
## F-statistic: 1.388e+04 on 1 and 1908 DF, p-value: < 2.2e-16datoscg$x05<-1/datoscg$age.yr**0.5 #cuando el modelo no corre por lo tanto busca otro valor
ajuste<-lm(y~x05,datoscg)
summary(ajuste)##
## Call:
## lm(formula = y ~ x05, data = datoscg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.42276 -0.14371 0.04043 0.18270 2.24037
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.23875 0.01915 221.3 <2e-16 ***
## x05 -6.16431 0.05887 -104.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2894 on 1908 degrees of freedom
## Multiple R-squared: 0.8518, Adjusted R-squared: 0.8517
## F-statistic: 1.096e+04 on 1 and 1908 DF, p-value: < 2.2e-16calcular la asintota
meter formulas
exp(4.23) #el valor de a viene de cuadro de ajuste anterior## [1] 68.71723ajuste
ajuste.nl<-nls (height.m~exp(a+b/(age.yr**c)),datoscg,start=c(a=4.3,b=-6,c=0.5))
summary(ajuste.nl)##
## Formula: height.m ~ exp(a + b/(age.yr^c))
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 4.18734 0.08937 46.85 <2e-16 ***
## b -7.61138 0.31873 -23.88 <2e-16 ***
## c 0.60362 0.03766 16.03 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.781 on 1907 degrees of freedom
##
## Number of iterations to convergence: 5
## Achieved convergence tolerance: 7.525e-06con esto nuestra curva quedaria de la siguiente forma
Ho(m)=e(4.18734-7.61138/t0.60362(años))
tener en cuenta la separacion de las curvas con el RMSE, ejemplo si se tiene una curva en 28,20 y otra 27, donde se tenga un RMSE DE 2.781 no se podra diferenciar ya que el rango de error es de 2.781, ´por lo qque ambas se sobreponen tendremos3 familias una por cada parametro
plot(datoscg$age.yr,datoscg$height.m) #FAMILIA a
meter formulas
despejar a
8=ea-7.6114/200.6036
log8=a-7.6114/20^0.6036
a=log8+7.6114/20^0,6036
a8<-log(8)+7.6114/20**0.6036
a13<-log(13)+7.6114/20**0.6036
a18<-log(18)+7.6114/20**0.6036
a23<-log(23)+7.6114/20**0.6036
a28<-log(28)+7.6114/20**0.6036
y8<-exp(a8-7.6114/x**0.6036)
y13<-exp(a13-7.6114/x**0.6036)
y18<-exp(a18-7.6114/x**0.6036)
y23<-exp(a23-7.6114/x**0.6036)
y28<-exp(a28-7.6114/x**0.6036)
plot(datoscg$age.yr,datoscg$height.m)
lines(x, y8, col="green")
lines(x, y13, col="green")
lines(x, y18, col="green")
lines(x, y23, col="green")
lines(x, y28, col="green")
esta familia es buena
FAMILIA b
meter formulas
despejar b
8=e4.1873-b/200.6036
log8=4.1873+b/20^0.6036
b=(log8-4.1873)*20^0,6036
b8<-(log(8)-4.1873)*20**0.6036
b13<-(log(13)-4.1873)*20**0.6036
b18<-(log(18)-4.1873)*20**0.6036
b23<-(log(23)-4.1873)*20**0.6036
b28<-(log(28)-4.1873)*20**0.6036
#para saber si vamos bien el b tiene que dar el valor de 7.61
y8<-exp(4.1873+b8/x**0.6036)
y13<-exp(4.1873+b13/x**0.6036)
y18<-exp(4.1873+b18/x**0.6036)
y23<-exp(4.1873+b23/x**0.6036)
y28<-exp(4.1873+b28/x**0.6036)
plot(datoscg$age.yr,datoscg$height.m)
lines(x, y8, col="red")
lines(x, y13, col="red")
lines(x, y18, col="red")
lines(x, y23, col="red")
lines(x, y28, col="red")
#FAMILIA C
meter formulas
despejar b
8=e4.1873-7.6114/20C
log8=4.1873-7.6114/20^c
7.6114/20^c=4.1416-log8
20^c=7.6114/(4.1873-log8)
clog8=log(7.6114/(4.1873-log8))
c=log(7.6114/(4.1873-log8))/log20
c8<-log(7.6114/(4.1873-log(8)))/log(20)
c13<-log(7.6114/(4.1873-log(13)))/log(20)
c18<-log(7.6114/(4.1873-log(18)))/log(20)
c23<-log(7.6114/(4.1873-log(23)))/log(20)
c28<-log(7.6114/(4.1873-log(28)))/log(20)
#para saber si vamos bien el C tiene que dar el valor de 0.60
y8<-exp(4.1873-7.6114/x**c8)
y13<-exp(4.1873-7.6114/x**c13)
y18<-exp(4.1873-7.6114/x**c18)
y23<-exp(4.1873-7.6114/x**c23)
y28<-exp(4.1873-7.6114/x**c28)
plot(datoscg$age.yr,datoscg$height.m)
lines(x, y8, col="blue")
lines(x, y13, col="blue")
lines(x, y18, col="blue")
lines(x, y23, col="blue")
lines(x, y28, col="blue")
#este familia es buena CURVAS ALGEBRAICAS
H1=a(1-e(-bt))c modelo original
se despajan H para dos mediciones
t1—-H1=a(1-e(-bt1))c t2—-H2=a(1-e(-bt1))c donde se obseva que lo unico que cambia es el tiempo, a,b y c son los mismos porque se trabaja con los mismos datos
despejar un parametro en t1 y sustituir en otro (t2)
a=H1/(1-e(-bt1))c se despejo a en H1
H2=H1(1-e(-bt2))c/(1-e(-bt1))c se sustituye en H2
El modelo da como resultado un conjunto de rectas con multiples asintotas
H2=H1(1-e(-bt2))c/(1-e(-bt1))c
despejar b en t1 anamorfica
H1/a=(1-e(-bt2))c
(H1/a)(1/c)=1-e(-bt1)
e(-bt1)=1-(H1/a)(1/c)
-bt1=log(1-(H1/a)^(1/c))
-b=log(1-(H1/a)^(1/c))/t1
Sustituir en el T2 anamorfica
H2=a(1-e(log(1-(H1/a)(1/c)))*(t2/t1))
H2=1-elog(1-(H1/a)(1/c)(t2/t1)c)
H2=a-(1-(H1/a)(1/c)(t2/t1)^c)
despeje de c anamorfica
H1/a=(1-e(-bt1)c)
log(H1/a)=c log(1-e^(-bt))
c=(log H1/a)/log(1-e^(-bt))
sustutuir
H2=a(1-e(-bt2)((log H1/a)/(log(1-e^(-bt1)))) se convierte en polimorfica con asintota unica
estos son tres ecuaciones ADAS