title: Testing the race model inequality in R

text: Lets say we have hypothetical participant that is presented with 3 stimuli conditions: a flashing dot ($C_x$), a sound blip ($C_y$) and combination of both ($C_z$). We ask this participant to respond as quickly as possible to the onset of the stimuli and we measure the reaction time (RT). It's been observed that responses are often observed to be substantially faster when both stimuli are presented simultaneously ($C_z$), compared to condition when they are presented separately ($C_x$ and $C_y$), and this effect has been labeled _redundant signals effect_ (e.g. Hershenson, 1962). Different models account for this effect, one of those models being _race model_ (Miller, 1982). Quoting Colonius et al (2006) definition: > Raab (1962) was the first to propose a race model for simple RT such that > (a) each individual stimulus elicits a detection process performed in > parallel to the others and (b) the winner’s time determines the observable > RT. This model suggests that RSE is generated by statistical facilitation: > If detection latencies are interpreted as (nonnegative) random variables, > the time to detect the first of several redundant signals is faster, on > average, than the detection time for any single signal. Urlich et al (2007) continues that: > According to race models, the observed RT distributions should satisfy, for > every value of $t$, the so-called _race model inequality_ (Miller, 1982): > > $F_z(t) \le F_x(t) + F_y(t), t > 0$ > > where $F_x$ and $F_y$ are the cumulative density functions (CDFs) of RT in > the two single-stimulus conditions $C_x$ and $C_y$, respectively, and $F_z$ > is the CDF of RT in the redundant-stimulus condition $C_z$. I wanted to test this model, but I found it hard to translate available Matlab code into R. Matlab code is available in the appendix of this paper by Urlich et al (2007), which can can be downloaded here, in case you don't have access to _Behavioural Research Methods_ journal. To summarize Urlich et al (2007) approach, his algorithm goes through the follwing steps: **Step 1: Estimating empirical cumulative density functions (CDFs) for every participant** (we got just one in the example here) and every stimulus condition ($C_x$, $C_y$, $C_z$). As Urlich et al (2007) points out: > these CDFs will be denoted as the _individual_ CDFs. Specifically, let > $G_x$, $G_y$, and $G_z$ be the individual estimates of two single-stimulus > CDFs, $F_x$ and $F_y$, and the redundant CDF, $F_z$, respectively. I thought I managed to compute $G_x$, $G_y$, and $G_z$, but I might be wrong:               #example reaction times for single hypothetical participants for 3 experimental  conditions     #this is the same dat set as one used by Urlich et al 2007     cx <- c(244, 249, 257, 260, 264, 268, 271, 274, 277, 291)     cy <- c(245, 246, 248, 250, 251, 252, 253, 254, 255, 259, 263, 265, 279, 282, 284, 319)     cz <- c(234, 238, 240, 240, 243, 243, 245, 251, 254, 256, 259, 270, 280)          #here I calculate quantiles corresponding to the given probabilities     #I use the default quantile algorithm (type 7 see R manual for 'quantile')     #it seems right, but I'm not 100% sure     psq <- seq(0,1,0.1)     gx <- quantile(cx, probs=psq, names=F)     gy <- quantile(cy, probs=psq, names=F)     gz <- quantile(cz, probs=psq, names=F)          #here I visualize the quantiles, looks more or less like in Urlich etal 2007     xyz <- data.frame(rt=c(gx,gy,gz),                        perc=rep(psq,3),                        cond=rep(c("gx(t)", "gy(t)","gz(t)"), each=length(psq)))     require(lattice)     xyplot(perc ~ rt, groups = cond, xyz, col=1, lty=3:1, type=c("g","l","p"),             key=list(lines=list(lty=3:1), text=list(labels=levels(xyz$cond))))      ![enter image description here](http://i.stack.imgur.com/tqIpc.png) **Step 2: From $G_x$ and $G_y$, one computes the bounding sum $B(t) = G_x(t) + G_y(t)$** (for each participant), which provides an estimate for each participant of the upper bound on the right side of _race model_ equation above. But then, I'm not sure how to compute the race model bounds, which seems trivial, but doesn't really work. The authors simply show the equation as , and I applied it in R as:               b = gx + gy      but then when I visualize it, it's wrong:               xyz <- data.frame(rt=c(gx,gy,gz,b),                        perc=rep(psq,4),                        cond=rep(c("gx(t)", "gy(t)","gz(t)","gx(t)+gy(t)"), each=length(psq)))     require(lattice)     xyplot(perc ~ rt, groups = cond, xyz, col=1, lty=4:1, type=c("g","l","p"),             key=list(lines=list(lty=4:1), text=list(labels=levels(xyz$cond))))      ![enter image description here](http://i.stack.imgur.com/scDbW.png) **Step 3 and 4: Percentile values are computed from $G_z$ and $B$** for each participant at certain prespecified probabilities that must be the same for all participants. Finally, the percentile values are aggregated across participants, which is omitted in Urlich etal (2007) example. I hope my explanation is a bit clearer now then it was before. **References** Miller, J. (1982). Divided attention: evidence for coactivation with redundant signals. Cognitive psychology, 14(2), 247–79. Ulrich, R., Miller, J., & Schröter, H. (2007). Testing the race model inequality: an algorithm and computer programs. Behavior research methods, 39(2), 291–302. Colonius, H., & Diederich, A. (2006). The race model inequality: interpreting a geometric measure of the amount of violation. Psychological review, 113(1), 148–54. doi:10.1037/0033-295X.113.1.148 Raab, D. H. (1962). Statistical facilitation of simple reaction times. Transactions of the New York Academy of Sciences, 24, 574-590. Hershenson, M. (1962). Reaction time as a measure of intersensory facilitation. Journal of Experimental Psychology, 63, 289-293.


title: Suggestions for multi-dimensional clustering

text: I am working in a genomics project and I ended up having a huge table with around 800 measurements (cases/rows), around 200 channels (columns/continuous variables) and 5 categories (one categorical column) I would like to do two things:   * Try to find sub-groups in the different levels of the categorical variable that I already have   * create a new classification of these 800 measurements based only in the information  I have been doing my homework and read about using different strategies like (k-means or PCA) but I have found that it is very useful to get rid of redundant variables. How can I choose these properly? Someone recommended me to use multinomial regression, any good resource you recommend to have a bite? I am using R. Many thanks


title: Nominal data multiple regression analysis

text: I have been given the following data and am supposed to build a multiple regression model to be able to predict the price of a car:   1. Price of Car - SPSS Measurement = Scale   2. Mileage of Car - SPSS Measurement = Scale   3. Age of Car - SPSS Measurement = Scale   4. **Number of Previous Owners - SPSS Measurement = Nominal**   5. Brand ID - SPSS Measurement = Nominal I am not sure if the Number of Previous Owners has been categorised properly, and if I should treat it as a dummy variable. Please advice! Thanks


title: Unit root for bounded variable

text: I have a time series which is discreet and bounded (between 0 and 100). It’s actually also increasing over time, from low to high. When I test for a unit root, I find it. And yet I wonder if unit root tests on bounded variables even make sense? And if unit-root testing doesn’t make sense in that case, what’s left of the whole cointegration framework? For instance, If I wanted to know if variable Y (0-100) is affected by another bounded variable X (say between 0-10), is it even possible to speak meaningfully about this?


title: Standard deviation when sample is significant part of population

text: I know that sample standard deviation and population standard deviation are different, and are used in different cases (first when N>>n, second when N=n). But what when N>n, although not by much? For example, n=7 and N=30. I guess none of these two would be appropriate, but is there another formula for such cases? Or at least a rule of thumb saying which of these two would be more correct?


title: How to compute marginals in Sum-Product Networks?

text: This should be fairly easy, but for some reason i'm having hard time getting it to work and I've spent a long time trying to figure it out myself. In the last paragraph of page 4 of the original Sum-Product Networks paper the authors described how to compute the posterior marginals of sum nodes, i.e $$P(Y_k = i | e) \propto w_{k,i} \frac{\partial S}{\partial S_K}\,,$$ where $Y_k$ is a sum node and $i$ is an edge to one of its child nodes. Let's assume that I've a very simple network like this one: ![Simple SPN](http://i.stack.imgur.com/Mxf1y.jpg) and I want to compute the marginal of the sum node given an evidence (which can be either: $e=\text{x}$ or $e=\text{not_x}$). The partial derivative of the sum node should be $1$ because it's the root: $$\frac{\partial S}{\partial S_{sum}} = 1$$ then the marginals should be: $$P(sum = 1 | e) \propto w_{sum,w_1} $$ $$P(sum = 2 | e) \propto w_{sum,w_2} $$ (where 1 and 2 referees to the edges) My problem is that these two equations don't depend on the evidence and they give the same results regardless of what the evidence is. I tried with different structures and I keep getting the same results. This is important because Sum nodes are seen as hidden variables in Sum- Product Networks and computing their posterior marginals is an important step in order to implement the Expectation-Maximization learning algorithm. How can I correctly compute the posterior marginal of the sum nodes?


title: If the likelihood principle clashes with frequentist probability then do we discard one of them?

text: In a comment recently posted here one commenter pointed to a blog by Larry Wasserman who points out (without any sources) that frequentist inference clashes with the likelihood principle. The likelihood principle simply says that experiments yielding similar likelihood functions should yield similar inference. Two parts to this question:   1. Which parts, flavour or school of frequentist inference specifically violate the likelihood principle?   2. If there is a clash, do we have to discard one or the other? If so, then which one? I will for the sake of discussion suggest that if we have to discard something then we should discard the parts of frequentist inference which clash, because Hacking and Royall have convinced me that the likelihood principle is axiomatic.


title: Obtaining the covariance matrix of a Gaussian random vector

text: I'm looking for some advice about a problem I've been assigned to solve. This is the problem. Suppose a car has to travel from P1 to P4, with intermediate points P2, P3. Say $X1$ is a r.v. that defines the travel time from P1 to P2 is normally distributed with $\mu$ = 25, $\sigma^2$ = 3, similarly $X2$ is a r.v. for the travel time from P2 to P3 and $X2 \sim N(15,3)$, finally X3 is a r,v, for the travel time from P3 to P4 and $X3 \sim N(20,1)$ . Then define a random vector $(T2,T3,T4)$ with the times when the car passes by P2,P3 and P4. I have to define the distribution of the random vector ${\bf t}=[T2,T3,T4]$ through its covariance matrix and its $\mu$ vector of expectations. Here is what I´ve done so far:   * $T2$ would be distributed just as $X1$. $T3$ would be defined as $X1+X2$, and $T4$ would be defined as $X1+X2+X3$. Then $T3$ would also be gaussian distributed since the sum of gaussians is a gaussian, and the distribution would be $N(25+15,3+3)=N(40,6)$. Following a similar reasoning, $T4 \sim N(25+15+20,3+3+1)=N(60,7)$. Then, if I define a vector ${\bf u}=[U_1,U_2,U_3]$ where $U_i \sim N(0,1)$, then ${\bf Au+\mu}$ should define the vector ${\bf t}$. If we define ${\bf A}$ as: $$ \begin{bmatrix} \sqrt3 & 0 & 0\\ 0 & \sqrt6 & 0\\ 0 & 0 & \sqrt7 \end{bmatrix} $$ and $\mu$ as $[15,40,60]$, then ${\bf Au+\mu}$ would define the distribution of vector ${\bf t}$, and the covariance matrix would be ${\bf \Sigma = AA^t}$. Is it correct? I am not sure because I think that somehow $T2$ would affect $T3$ and $T4$ and therefore $Cov(T2,T3), Cov(T3,T4)$ would not be zero. Also I have to find the probability that the travel time from $P1$ to $P4$ is greater than 70. I think this would be solved as $P(T4>70)$, since $T4 \sim N(60,7)$,if this is true I just simply have to normalize and take a look and the normal distribution table or go to R and type pnorm(70,60,7). Do you think it is correct? Finally, I have to estimate the moment when the car passed by P3, given that it is known that it arrived at P4 in 65.234 time units. Also, if I know that it passed by P2 in 26.1 time units, I have to estimate the moment when the car passed by P3. In this point I'm lost and don't know exactly what to do. So this is it, sorry for the looong post. Thanks for reading. Any comments and suggestions are deeply apreciated.


title: Standard deviation for values below 1

text: When analyzing a dataset based on percentages, I have on occasions worked with data as "full" value (i.e. "50") or "reduced" value (i.e. ".50"). However, it just occurred to me that this could have a serious impact on how standard deviation and variance. If my standard deviation and variance are above 1, the standard deviation will be smaller than the variance. But if they are below 1, the standard deviation will be **_bigger_** than the variance. Does it mean it is good practice to never have a dataset resulting in a variance below 1? Or is there a way to account for this?


title: What is cross section in "cross section of stock return"?

text: Can somebody give the definition of cross section in "cross section of stock return"? Thanks


title: Deciding on power to use for box-cox transformation

text: In R, I am using the `bc` function to do a box-cox transformation. What factors do I need to consider when setting `p` (the power argument)?


title: How to test a null hypothesis regarding 2 dichotomous variables?

text: After an interesting lecture on Bayesian statistics (research master), all 4 males went immediately to the (female) teacher to ask questions. The 8 females did not. I am interested in how to test a two-tailed significance null hypothesis with this, in `R`. H0: no difference between males and females in talking to teacher yes/no H1: difference between males and females in talking to teacher yes/no What is the correct test? I have a difficulty with this since:   * unequal group length   * should we use 2 variables? (male/female and outcome yes or no), or just one which combines? I tried many things, but I can't decide what is the correct way. ( I understand that this is kind of a ridiculous question, and understand the irony for such a test after having had Bayesian statistics, but I am really interested)


title: How to visualise a 3D data-set for someone over the Internet?

text: I was just about to ask a question concerning a certain 3D data-set (each point in space has a value associated), but couldn't figure out how to visualize it for you. I can understand it myself by rotating the point set in a Matlab figure, but it would be difficult to grasp from just seeing a screenshot: ![ ](http://i.stack.imgur.com/nPzeD.png) (marker size and color are both proportional to the value at that point in space) What's the best way of solving this problem?


title: Variation on binomial experiment

text: Consider drawing stones from an urn with replacement. Suppose the urn contains a large number of stones, and each stone has a color, and the set of possible colors is large. Suppose $N$ stones are drawn one at a time from the urn, and replaced each time. Suppose $S$ of them have the same color, $C$. Can I write a distribution telling what I know about the fraction $F$ of stones having color $C$ in the urn? Clearly it depends on my prior knowledge. Suppose I "know nothing" - i.e. all colors are equally likely. I'm tempted to say $F \sim \beta(S+1, (N-S)+1)$ as if I knew in advance what color I were looking for, but I didn't. Should I instead discount one of the $C$-colored samples and treat it as having given me prior knowledge, decreasing both $S$ and $N$ by $1$ ?


title: Combination of 2 random variables that is perfectly correlated with another random variable

text: I'd appreciate some pointers for the following question: Given 3 random variables $X_i$, $i = 1, 2, 3$ such that none of them have a pairwise correlation of $+1$ or $−1$. Can we find a combination of $X_1$ and $X_2$ that has a correlation of $+1$ with $X_3$?


title: Heterogeneity with two studies

text: With data from two centres I want to account for potential heterogeneity or confounders between two centers. So the analysis will initially be stratified by clinical center and a chi square test performed with one degree of freedom. Is this appropriate with just two centres? Or is there an alternative?


title: Combining LASSO coefficients across imputed datasets

text: I am using the LASSO with multiple imputed datasets and I am not sure how should I combine the coefficients obtained on the different imputed datasets. I could simply average them (as I would do had I computed the coefficients using ordinary least squares), but then since the set of variables with non- zero coefficients are not necessarily equal on all imputed datasets a variable would end up in my model even if it only appears on one of the imputed datasets. Another idea would be to take the average but only if the coefficient is different from zero in at least half (or some other proportion) of the imputed datasets. I would be really thankful any suggestions or references on how to do this (I do not know if this matters but I am selecting the LASSO tuning parameter on each imputed dataset using cross-validation).


title: Specifying parameter constraints in nls()

text: Is it possible to specify that one parameter must be larger than another parameter in an nls call in my R script? Here's my nls call:               fit <- nls(y ~ ifelse(g, m1 * (x - x0) + y0, m2 * (x - x0) + y0),                 start = c(m1 = -1, m2 = 1, y0 = 0, x0 = split),                 algorithm = "port",                 lower = c(m1 = -Inf, m2 = -Inf, y0 = -Inf, x0 = split),                 upper = c(m1 = Inf, m2 = Inf, y0 = Inf, x0 = (split+1)),                 data=data.frame(x,y))      It is finding values for the parameters m1, m2, x0, and y0. But I want to require that m2 must be LARGER than m1. How can I do this?


title: How to illustrate the "Coefficient of determination" using graphics (preferably using R)?

text: I came across this nice illustration of the "Coefficient of determination" (source): ![enter image description here](http://i.stack.imgur.com/ukljW.png) Which leads me to ask two questions:   1. How to do it with R? (I guess my main question would be how to deal with the opacity)   2. Are there other interesting/useful visualizations of the "Coefficient of determination"? (and again, is there an easy way to make them in R?) Thanks.


title: R time-series forecasting with neural network, auto.arima and ets

text: I've heard a bit about using neural networks to forecast time series. How can I compare, which method for forecasting my time-series (daily retail data) is better: auto.arima(x), ets(x) or nnetar(x). I can compare auto.arima with ets by AIC or BIC. But how I can compare them with neural networks? For example:                  > dput(x)      c(1774, 1706, 1288, 1276, 2350, 1821, 1712, 1654, 1680, 1451,       1275, 2140, 1747, 1749, 1770, 1797, 1485, 1299, 2330, 1822, 1627,       1847, 1797, 1452, 1328, 2363, 1998, 1864, 2088, 2084, 594, 884,       1968, 1858, 1640, 1823, 1938, 1490, 1312, 2312, 1937, 1617, 1643,       1468, 1381, 1276, 2228, 1756, 1465, 1716, 1601, 1340, 1192, 2231,       1768, 1623, 1444, 1575, 1375, 1267, 2475, 1630, 1505, 1810, 1601,       1123, 1324, 2245, 1844, 1613, 1710, 1546, 1290, 1366, 2427, 1783,       1588, 1505, 1398, 1226, 1321, 2299, 1047, 1735, 1633, 1508, 1323,       1317, 2323, 1826, 1615, 1750, 1572, 1273, 1365, 2373, 2074, 1809,       1889, 1521, 1314, 1512, 2462, 1836, 1750, 1808, 1585, 1387, 1428,       2176, 1732, 1752, 1665, 1425, 1028, 1194, 2159, 1840, 1684, 1711,       1653, 1360, 1422, 2328, 1798, 1723, 1827, 1499, 1289, 1476, 2219,       1824, 1606, 1627, 1459, 1324, 1354, 2150, 1728, 1743, 1697, 1511,       1285, 1426, 2076, 1792, 1519, 1478, 1191, 1122, 1241, 2105, 1818,       1599, 1663, 1319, 1219, 1452, 2091, 1771, 1710, 2000, 1518, 1479,       1586, 1848, 2113, 1648, 1542, 1220, 1299, 1452, 2290, 1944, 1701,       1709, 1462, 1312, 1365, 2326, 1971, 1709, 1700, 1687, 1493, 1523,       2382, 1938, 1658, 1713, 1525, 1413, 1363, 2349, 1923, 1726, 1862,       1686, 1534, 1280, 2233, 1733, 1520, 1537, 1569, 1367, 1129, 2024,       1645, 1510, 1469, 1533, 1281, 1212, 2099, 1769, 1684, 1842, 1654,       1369, 1353, 2415, 1948, 1841, 1928, 1790, 1547, 1465, 2260, 1895,       1700, 1838, 1614, 1528, 1268, 2192, 1705, 1494, 1697, 1588, 1324,       1193, 2049, 1672, 1801, 1487, 1319, 1289, 1302, 2316, 1945, 1771,       2027, 2053, 1639, 1372, 2198, 1692, 1546, 1809, 1787, 1360, 1182,       2157, 1690, 1494, 1731, 1633, 1299, 1291, 2164, 1667, 1535, 1822,       1813, 1510, 1396, 2308, 2110, 2128, 2316, 2249, 1789, 1886, 2463,       2257, 2212, 2608, 2284, 2034, 1996, 2686, 2459, 2340, 2383, 2507,       2304, 2740, 1869, 654, 1068, 1720, 1904, 1666, 1877, 2100, 504,       1482, 1686, 1707, 1306, 1417, 2135, 1787, 1675, 1934, 1931, 1456)      Using auto.arima:               y=auto.arima(x)     plot(forecast(y,h=30))     points(1:length(x),fitted(y),type="l",col="green")      ![enter image description here](http://i.stack.imgur.com/uwSqY.png)               > summary(y)     Series: x      ARIMA(5,1,5)                              Coefficients:              ar1      ar2     ar3      ar4      ar5      ma1     ma2      ma3     ma4      ma5           0.2560  -1.0056  0.0716  -0.5516  -0.4822  -0.9584  1.2627  -1.0745  0.8545  -0.2819     s.e.  0.1014   0.0778  0.1296   0.0859   0.0844   0.1184  0.1322   0.1289  0.1388   0.0903          sigma^2 estimated as 58026:  log likelihood=-2191.97     AIC=4405.95   AICc=4406.81   BIC=4447.3          Training set error measures:                        ME     RMSE      MAE       MPE     MAPE      MASE     Training set 1.457729 240.5059 173.9242 -2.312207 11.62531 0.6157512      Using ets:               fit <- ets(x)     plot(forecast(fit,h=30))     points(1:length(x),fitted(fit),type="l",col="red")      ![enter image description here](http://i.stack.imgur.com/9UngX.png)                > summary(fit)      ETS(M,N,N)            Call:       ets(y = x)              Smoothing parameters:          alpha = 0.0449              Initial states:          l = 1689.128              sigma:  0.2094                AIC     AICc      BIC       5570.373 5570.411 5577.897            Training set error measures:                         ME     RMSE      MAE      MPE     MAPE      MASE      Training set 7.842061 359.3611 276.4327 -4.81967 17.98136 0.9786665      In this case auto.arima fits better then ets. Let's try sing neural network:                library(caret)      fit <- nnetar(x)      plot(forecast(fit,h=60))      points(1:length(x),fitted(fit),type="l",col="green")      ![enter image description here](http://i.stack.imgur.com/M8HIT.png) From the graph, I can see, that neural network model fits quite well, but how can I compare it with auto.arima/ets? How can I compute AIC? Another question is, how to add confidence interval for neural network,if it is possible, like it is added automatically for auto.arima/ets.? Any help and advises would be really appreciated.


title: What are the differences between document classification and clustering when working with a single topic?

text: I am doing some web page clustering work and I'm going to use cosine similarity as my distance measure. Even though cosine similarity is a clustering technique, I have to give training data in order to build the query vector. Clustering algorithm doesn't need training data in the sense of with labeled classes, but how do you build the query vector if you don't give the training data in the cosine similarity calculation? I am only interested in a single topic (sports) so if I do it with 2 clusters, when a new document is fed, if it is clustered to the cluster 1 (say sports), then I'll take that document or else it will be rejected. In this case how this differs from single-class classification?


title: How to assess regression performance when using PCA components as response variables?

text: I have a regression problem for multi-dimensional response, say like 1000 dimension response, without using PCA to the response, my mean square error for regression is like ~80. Then, I apply PCA to these 1000 dimensions, and use 3 of them which represent 90% variation of the response, and perform the regression again. Now the mse is ~8. So how can I interpret this results? Can I claim that using 3 PCA components can get better regression since mse is 8 which is significant smaller than 80?


title: Post stratification weights in survey package in R

text: I have data collected from a survey administered on a subset of the population. I also have the population proportions of variables such as gender, race and housing type. I would like to combine the weights from each separate cross tab (of gender, race and housing type) such that the weighted proportions of my survey data matches that of the population. I have tried the following:               library(survey)          gender.population <- read.table("http://dl.dropbox.com/u/822467/Gender.csv", header = TRUE, sep = ",")     housing.population <- read.table("http://dl.dropbox.com/u/822467/Housing.csv", header = TRUE, sep = ",")     race.population <- read.table("http://dl.dropbox.com/u/822467/Race.csv", header = TRUE, sep = ",")     survey.sample <- read.table("http://dl.dropbox.com/u/822467/survey.sample.csv", header = TRUE, sep = ",")          survey.object.sample <- svydesign(id = ~1, data = survey.sample)     survey.object.sample.weighted <- rake(survey.object.sample, list(~gender, ~housing, ~race), list(gender.population, housing.population, race.population))          str(survey.object.sample.weighted$postStrata)      I see from survey.object.sample.weighted$postStrata that weights have been assigned separately for each of the variable. My question is: Is it possible to get 1 weight for each subject instead of 3 weights as shown in the package?


title: Arima modelling estimate in SAS and find out if model is adequate or not

text: I have been attempting to do an ARIMA modelling in SAS. The series is not stationary but when I estimate the ACF and PACF p-values I don't get appropriate answers to find out if my model is adequate or not.Can anyone help me as may be I have been choosing the wrong p and q values. year age-group members.For stationarity I have used Sas which shows me that series is stationary with constant mean and variance.![ model showing Stationary series for members](http://i.stack.imgur.com/GwpIk.png) year age-group members 1997 16to21 21.1 1997 22to29 47.5 1997 30to39 61.7 1997 40to49 65.1 1997 50to54 64.1 1997 55to59 58.7 1997 60to64 48.2 1997 65ANDOVER 12.2 1998 16to21 18.6 1998 22to29 46.3 1998 30to39 61.1 1998 40to49 65.4 1998 50to54 64 1998 55to59 59 1998 60to64 45.4 1998 65ANDOVER 12.5 1999 16to21 19.6 1999 22to29 47.3 1999 30to39 61.6 1999 40to49 65.8 1999 50to54 64.8 1999 55to59 58.5 1999 60to64 45.2 1999 65ANDOVER 10.1 2000 16to21 20.9 2000 22to29 47.6 2000 30to39 61.5 2000 40to49 65.8 2000 50to54 64.6 2000 55to59 58.8 2000 60to64 44 2000 65ANDOVER 10.7 2001 16to21 20.7 2001 22to29 47.1 2001 30to39 61.7 2001 40to49 65.7 2001 50to54 64.5 2001 55to59 59.2 2001 60to64 42.6 2001 65ANDOVER 11.7 2002 16to21 19.4 2002 22to29 49.2 2002 30to39 62.8 2002 40to49 67.1 2002 50to54 65.8 2002 55to59 61.2 2002 60to64 45.2 2002 65ANDOVER 10.1 2003 16to21 17.3 2003 22to29 45.2 2003 30to39 59.7 2003 40to49 65.1 2003 50to54 64.9 2003 55to59 59.1 2003 60to64 44.5 2003 65ANDOVER 10.7 2004 16to21 14.1 2004 22to29 44 2004 30to39 60.5 2004 40to49 66.2 2004 50to54 66.3 2004 55to59 61.4 2004 60to64 45.9 2004 65ANDOVER 12.1 2005 16to21 11.4 2005 22to29 38.7 2005 30to39 56.1 2005 40to49 62.3 2005 50to54 63.2 2005 55to59 58 2005 60to64 42.7 2005 65ANDOVER 6.8 2006 16to21 11.7 2006 22to29 38.1 2006 30to39 56 2006 40to49 62.1 2006 50to54 62.8 2006 55to59 58.8 2006 60to64 44.7 2006 65ANDOVER 10 2007 16to21 11.4 2007 22to29 36.9 2007 30to39 54.8 2007 40to49 60.7 2007 50to54 62.3 2007 55to59 58.7 2007 60to64 45.1 2007 65ANDOVER 10.9 2008 16to21 9.2 2008 22to29 34.8 2008 30to39 53.6 2008 40to49 60 2008 50to54 62.9 2008 55to59 59 2008 60to64 45.4 2008 65ANDOVER 12.3 2009 16to21 9.6 2009 22to29 35.4 2009 30to39 53.2 2009 40to49 60.3 2009 50to54 62.4 2009 55to59 59.7 2009 60to64 46.3 2009 65ANDOVER 12.5 2010 16to21 9.1 2010 22to29 34.1 2010 30to39 52.6 2010 40to49 59.6 2010 50to54 62 2010 55to59 59.2 2010 60to64 46.3 2010 65ANDOVER 13.9 2011 16to21 7.7 2011 22to29 31.8 2011 30to39 50.7 2011 40to49 57.9 2011 50to54 60.1 2011 55to59 57.3 2011 60to64 45 2011 65ANDOVER 14.6 2012 16to21 6.5 2012 22to29 30.2 2012 30to39 49.1 2012 40to49 56.1 2012 50to54 58.2 2012 55to59 55.5 2012 60to64 43.9 2012 65ANDOVER 14.8


title: Monte Carlo simulation exercise

text: The assignment is to run a simulation on excel to figure out the probabilities on obtaining a certain goal: John and Jane Doe are planning to save money to pay a house for their 6-month- old son, Patrick. They have decided that they would like to have $500,000 saved by the time Patrick is ready for college 17 years from today: John and Jane are planning to save $20,000 (at the beginning of) each year (ignore taxes). Assume that the return each year varies and is sampled from the same normal distribution (mean rate of return is 4%, standard deviation is also 10%, i.e. one standard deviation away is either a -6% return or an 14% return). (Hence in each trial there are 17 different rates of return, one for each year.) Run a simulation of this investment strategy with at least 1000 trials. Here's my attempt on the problem (click on the photo to reveal a larger version): ![Excel Spreadsheet of Simulation](http://i.stack.imgur.com/0ewt8.png) I am not to sure on how to continue from there. Am I missing other components to solve the problem? Is 'n' the number of trials? Any advice would be appreciated!


title: Explaining quantile regression to nonstatisticians

text: I recently submitted a paper, in which I used quantile regression, to a psychology journal. Although I thought I had already put enough thought in a clear exposition of quantile regression, the reviewers asked for better explanations of the quantile regression technique being only familiar with standard OLS regression. So, what is the best way to explain quantile regression, in an empirical paper, to non-statisticians?


title: How many features to overfit the classification?

text: Recently, I have got some 'strange' comments from the reviewer of my paper. In my paper, I discussed a novel feature extraction method, and then I compared three classification methods for my binary classification task. I applied kNN, Neural Networks, and Support Vector Machine classifiers. Simply speaking I have got: 2 Class labels (Binary classification labelled with either 1 or 0); **50** Cases (20 labelled as 1 and 30 labelled as 0); 18 Features extracted (some physical quantities); **12** Features remained after feature selection; After the feature selection, I have 12 features remained and the classifiers were running on the features with a leave-one-case-out cross validation. For the results I have got AUC of kNN: 0.91 Neural Networks: 0.92 SVM: 0.97 My conclusion is that SVM outperformed other classifiers (and also I have got other quantitative results). And the overall performance proved that my feature extraction method is robust and suitable for my binary classification problem. The reviewer's comments are:   1. The authors applied three different learning strategies for XXXX classification. Even if these were all equally accurate one must outperform the others **purely by chance**.   2. 12 features used for 50 test cases. Will the features **overfit** the classifier? a. For the first question, I think the reviewer is correct, the SVM might outperform the other two by chance. However, isn't it normal in comparison of classifiers? The classification performance is data dependent. For example, in real application, if there is an outlier, the classification performance can be degenerated. Also, kNN can outperform SVM in some case. Or for another example, if we take out the cases kNN mis-classified. We may get equal performance of the kNN and SVM. In summary, I think one classifier outperform others can be 'by chance'. And the the experimental accuracy is likely to be greater than the true accuracy. What do you guys think? b. For the second question, I don't understand. Is there a relationship between the number of features and overfitting? In an extreme case: can I extract 100 features to classify 50 cases? Any theory behind for the number of features? (please list references for your thought if you can find some.) Thanks very much. A


title: How robust is ANOVA when group sizes are unequal and residuals are not normally distributed?

text: I understand that there are degrees of trust one can have in the output of ANOVA. However, naturally I want to maximize the amount I can trust my results. **Question:** If I have data that violate the assumption of normality (of residuals), how robust is a factorial ANOVA when cell sizes are unequal? My cell sizes range from 42 to 55, for a 2x3 between-subjects ANOVA.


title: Whitening Transformation using a Hadamard product Variance Matrix

text: I want to whiten a vector $X$ by transforming the variance-covariance matrix so the variance-covariance matrix of the transformed series will be the identity matrix $I$. $X$ is a time-series column vector with variance-covariance matrix $M$ and mean $0$. If $X$ was a set of correlated data, then the covariance matrix $M$ could be written as the outer product of $X$ with itself: $M=\mathrm{E}\left[XX^T\right]$. $M$ is thus a variance-covariance matrix, symmetric and positive semidefinite, has a square root $M^{{1}/{2}}$, so $M^{{1}/{2}}\left(M^{{1}/{2}}\right)^T=M$. If $M$ is positive definite then $M^{{1}/{2}}$ is invertible, so the vector $Y=M^{-{1}/{2}}X$ has covariance matrix: $\mathrm{Cov}\left[Y\right]=\mathrm{E}\left[YY^T\right] = M^{-{1}/{2}}\mathrm{E}\left[XX^T\right]\left( M^{-{1}/{2}}\right)^T=M^{-{1}/{2}}M\left( M^{-{1}/{2}}\right)^T=\left(M^{-{1}/{2}}M^{{1}/{2}}\right)\left(M^{{1}/{2}}\right)^T\left( M^{-{1}/{2}}\right)^T=I$ However, when $X$ is a time series $M$, the variance-covariance matrix is constructed from the autocovariance function of $X$ as: $\left(\begin{array}{cccc} \gamma\left(0\right) & \gamma\left(1\right) & \cdots & \gamma\left(n-1\right)\\ \gamma\left(1\right) & \gamma\left(0\right) & \cdots & \gamma\left(n-2\right)\\ \vdots & \vdots & \ddots & \vdots \\ \gamma\left(n-1\right) & \gamma\left(n-2\right) & \cdots & \gamma\left(0\right) \end{array}\right)$ which is not the same as $M=\mathrm{E}\left[XX^T\right]$. However, this can be written as the Hadamard product $W \circ M$, where $W$ is a known set of weights. When I apply the method, I can find a square root of $W \circ M$ that decorrelates the dataset $X$. However, the decorrelated dataset does not have unit variance - the varaince is lower because of the effect of $W$. However, after doing this for around 100 different datasets, the required correction factor is 0.75, or thereabouts. Is there a way to derive this correction factor directly from the matrix algebra? Many thanks in advance. Just for the record, this question is not being asked in order to answer to any test or assignment.


title: Confounder - definition

text: According to M. Katz in his book Multivariable analysis (Section 1.2, page 6), " _A confounder is associated with the risk factor and causally related to the outcome._ " Why must the confounder be _causally_ related to the outcome? Would it be enough for the confounder to be _associated_ with the outcome?


title: Is $H=\min(t_1,...,t_n)$ a Copula?

text: Please help me prove the following: n-Box is defined as $B=[a_1,b_1]\times[a_2,b_2]\times[a_3,b_3]\times...\times[a_n,b_n]$ Cartesian product of $n$ closed intervals, where $a_i$ and $b_i$ are all from $R$ and $a_i\le b_i$ for $i=1,...,n$. Vertices of the box are $c=(c_1,...,c_n)$ where $c_k=a_k$ or $c_k=b_k$. We have a function defined as $H=\min(t_1,...,t_n)$. We define $H$ volume of the box as $V_H(B)=\sum \mathrm{sgn}(c)H(c)$ where the sum is taken over all vertices of the box. $\mathrm{sgn}(c)=1$ if $c_k=a_k$ for an even number of $k$'s $\mathrm{sgn}(c)=-1$ if $c_k=a_k$ for an odd number of $k$'s Prove that $V_H(B)=\max(\min(b_1,...,b_n)-\max(a_1,...,a_n),0)$ That would also be a proof that $H$ is Copula. Thanks in advance.


title: Having a conjugate prior: Deep property or mathematical accident?

text: Some distributions have conjugate priors and some do not. Is this distinction just an accident? That is, you do the math, and it works out one way or the other, but it does not really tell you anything important about the distribution except for the fact itself? Or does the presence or absence of a conjugate prior reflect some deeper property of a distribution? Do distributions with conjugate priors share some other interesting property or properties that other distributions lack that causes those distributions, and not the others, to have a conjugate prior?


title: Fisher's method for combing p-values - what about the lower tail?

text: I have a bunch of independent p-values and now I want to combine them using the Fisher's method. Each of the individual p-values is coming from a one- sided test. I am just a little bit confused about the "side" of the Fisher's method test, i.e. when I calculate the Fisher's method p-value in R, I use:               1 - pchisq( -2*sum(log(p-values)), df)      where `df = 2*length(p-values)`. Is this a one-sided test? It should be because when the test statistic `-2*sum(log(p-values))` is much smaller than `df`, then the Fisher's p-value is close to 1. There should be a problem here, right? How shall I (or should I?) reject the null if my test statistic is very small? I am just uncomfortable with the "close-to-1" p-values. BTW, I use the method for testing model goodness-of-fit, and a small p-value is to indicate model lack-of-fit. Thanks!


title: Statistical Significance test for difference in FScore

text: For a classification task, I have developed two methods and FScore(harmonic mean of precision and recall) of both the classes serves as the performance criteria. How can I check whether the difference between the performance of the two systems is statistically significant or not. For eg. I was reading a paper by Navigli and in the evaluation he cites that "all the differences between 2 systems are statistically significant (p < 0.01)" I would like to know how to perform this test of statistical significance. EDIT : No clue is there regarding which test is been talked about. Is there some standard test used in NLP community? (the work is related to Natural Language Processing)


title: Is highly correlated factors in a prediction model a problem?

text: I have built a logistic regression model with two or more highly correlated factors. I did this by doing a bagging procedure. In my understanding having highly correlated factors in a _prediction_ model is not an issue (especially after bagging). But in a model where you are trying to _explain_ the target, this could be an issue. I remember something like the above was stated in a machine book but I can't find it. Do you guys agree with this viewpoint and can point me to a reference or two on this?


title: Repeated measures ANOVA in R for non-orthogonal design

text: I have data from a two-factor within-subjects experiment design where the conditions are not orthogonal. Factor one (Location) has three levels; factor two (Stimulus) has three levels, one of which is "no stimulus". When multiplying these factors I get nine conditions, three of which are equivalent. All three Location x "no stimulus" conditions are equivalent because of lack of stimulus. As a result I only have observations from 7 conditions (just a single Location x "no stimulus" condition, intended to be used as a control condition). I just chose an arbitrary level of Location for that condition for the sake of coding the data. My question is - how can I analyse my data in R using ANOVA? I've tried building a model with the lme() function in the nlme package (following a textbook example for orthogonal designs), but I get an error when trying to build the model (presumably because of missing conditions?):               > model <- lme(Y ~ Location * Design, data, random = ~ 1 | Subject / Location / Design, method="ML")     Error in MEEM(object, conLin, control$niterEM) :        Singularity in backsolve at level 0, block 1      The lmer() function in the lme4 package also gives an error:               model <- lmer(Y ~ Location * Design + (1|Subject) + (1|Location) + (1|Design) + (1|Location:Design), data, REML=FALSE)     Error in mer_finalize(ans) : Downdated X'X is not positive definite, 8.      The ezANOVA() function in the ez package also gives an error. Any advice on how I could approach this analysis would be greatly appreciated! Would it be bad to duplicate the "no stimulus" observations for each of the two Location x "no stimulus" conditions which don't have observations?


title: Cross-Correlating more than Two Time Series

text: this is my first post to stackexchange, I have read many solutions, especially to programming questions in R and have appreciated their helpfulness. Currently I find myself stuck with a time series analysis problem. I am an aspiring dendrochronologist, measuring the annual tree-rings from a particularly difficult tropical tree. I have $n=37$ stem cross-sections, that I will soon finish measuring the area of, and form time series that estimate tree growth per year. Each of these stems have between 60 and 80 years in their series, but have unknown dates for each tree-ring (many died in unknown years). I need to cross-correlate all the tree samples together and adjust each series to the best lags, before averaging all the series into an estimated annual forest growth for the population. I would like to trend away from traditional methodology, of plotting each time series on a graph and visually cross-dating, and prefer to use a robust function in R, python, c++, or matlab. Currently I have written two algorithms in R, that both rely on the cross- correlation function _ccf() acf{stats} package_ which can return the pairwise euclidean distance between all samples. These pairwise correlations are used for forming correlation matrices between all samples, with a predicted lag for each tree growth series. When looking at a high $\text{lag.max}=40$, you cannot test every possible permutation of even $n=20$ samples because this becomes: $$(2*\text{max.lag}+1)^{(n-1)} = \text{ states; or }(2*40+1)^{(20-1)}=1.8\text{e}36 \text{ states}$$ which is not testable in a reasonable amount of time. I circumvented testing all possibilities by binning samples into groups of four, and solving all possibilities within these bins, then cross-correlating these bins to each other and adjust the series in the bins accordingly. This becomes slow as you add more samples, even though its only testing a small subset of all possibilities, which potentially sacrifices the optimal solution. **I am wondering** if there is a good algorithm that can find the optimal cross-correlation between all time series that can avoid testing all 1.8e36 possibilities. I have made some attempts to use machine learning algorithms inside of Weka, but am unsure how they could be easily applied to this problem. Thank you for reading, Brook


title: Why no mention of penalized regression techniques in Applied Logistic Regression, 3rd edition, by Hosmer, Lemeshow, and Sturdivant?

text: Just ordered this textbook, and Wow, the complete omission of this subject from an otherwise excellent reference on logistic regression is a bit surprising. The 2nd edition was published in 2000 - there's been a lot of research into penalized regression techniques since then: for example, the lasso and elastic net. I know this text is aimed at epidemiologists and public health specialists. Have researchers in these fields ignored any developments that smack of data mining?


title: How to compare observed vs. expected events?

text: Suppose I have one sample of frequencies of 4 possible events:               Event1 - 5     E2 - 1     E3 - 0     E4 - 12      and I have the expected probabilities of my events to occur:               p1 - 0.2     p2 - 0.1     p3 - 0.1     p4 - 0.6      With the sum of the observed frequencies of my four events (18) I can calculate the expected frequencies of the events right?               expectedE1 - 18 * 0.2 = 3.6     expectedE2 - 18 * 0.1 = 1.8     expectedE1 - 18 * 0.1 = 1.8     expectedE1 - 18 * 0.6 = 10.8      How can I compare observed values vs expected values? to test if my calculated probabilities are good predictors? I thought of a chi-squared test, but the result change with the sample size (n=18), I mean, if I multiply observed values by 1342 and use the same method the result is different. Maybe a wilcox paired test works, but what do you suggest? If can suggest in R, it would be better.


title: Can someone explain this confidence interval question to me?

text: What proportion of shoppers at a large appliance store actually make a large ticket purchase? To estimate this proportion within 10% and be 95% confident of the results, how large a sample should you take? Assume that you have no idea what proportion of all shoppers actually make a large ticket purchase A) 68 B) 97 C) 49 D) 83


title: Two uneven groups (101 vs 21) measured twice; what repeated measures ANOVA equivalent can help?

text: I have a question that I have difficulties in answering: I have two groups of patients tested twice at the same continuous variable. Group A=104 & GroupB=21. So two groups, one continuous variable, measured twice. Due to the uneveness of the samples I think that repeated measures ANOVA is not suitable. Is there another way of testing the patient groups differences? Maybe a non- parametric equivalent? I prefer an SPSS solution, but I could use R as well. Thanks in advance.


title: How to create a bivariate normal distribution

text: I have a data set like this:               df          Income     Education_in_years     40,000     10     50,000      9     70,000     12     30,000      5     100,000    20      I would like to create a bivariate distribution from this and try to guess probability of income given eduction in years. I can build the linear model as follows:               lin <- lm(Income~Eduction_in_years, data=df)      I could come up with a formula like this:               Income = a*education_in_years + e      What I would like to do is, given the model, create a bivariate normal, and run a 'what-if' analysis to determine `income` given `eduction_in_years`. Can somebody walk me through how to create a bivariate normal distribution from this dataset and determine the income level given the education?


title: linear (?) mixed models help?! - simple question from a simple mind

text: I have a stunningly simple data set and, unfortunately, an equally simple mind when it comes to stats. This is how my data are set up: I asked people to sit small exams on different days of the week (Friday, Saturday, Sunday and Monday) and their performances in these tests have been calculated as a percentage. Some of them only sat the test on one day, others on more than one day and as such I have a data set that can neither be analysed with a paired nor independent samples test without cutting out data (and my sample size is small enough so that's a no go). A mixed model is what I have been told I need and as my data are normally distributed am I correct in saying that a linear mixed model is the one to choose? My study question is whether score is influenced by day of the week so I understand that I set up the model with 'score' as dependent variable, 'ID' as random effect and 'day' as fixed effect - please tell me if that is wrong! So [1] can anyone point me in the direction of where to get a walkthrough for a LMM (or alternative test if that is incorrect), [2] tell me how to calculate the effect sizes, and [3] what exactly do you report from the wall of tables and figures that spews out of SPSS when it has completed the calculations?! All help much appreciated! * * * PS - I have just attempted to tag this post and discovered there is a different between mixed-models and mixed-effects-models. I did not know that so if anyone can explain that also, that would also be swell.


title: Step-by-step example of predicting time series with ARIMAX or ARMAX model?

text: Could someone give me a step-by-step example of time series prediction using ARIMAX or ARMAX model? The example doesn't need to be long or complicated. It could be for example forecasting temperature with past data of only ten values (e.g [15, 16, 17, 15, ..., 12]) You don't have to calculate autocorrelation etc. step-by-step but I need it to be clear A) What are you doing and why, B) Which data you're using to get what etc.. For example, I would like to see the procedure of finding AR-, MA-, and X- lags, possible integrations, finding the error terms $y_i - \widehat{y}_i = e_i$, solving the coefficients for AR-, MA- and X-parts and making the predictions.


title: How does the correlation coefficient differ from regression slope?

text: I would have expected the correlation coefficient to be the same as a regression slope (beta), however having just compared the two, they are different. How do they differ - what different information do they give?


title: Can Random Forests do much better than the 2.8% test error on MNIST?

text: I haven't found any literature on the application of Random Forests to MNIST, CIFAR, STL-10, etc. so I thought I'd try them with the _permutation-invariant_ MNIST myself. In **R** , I tried:               randomForest(train$x, factor(train$y), test$x, factor(test$y), ntree=500)      This ran for 2 hours and got a 2.8% test error. I also tried **scikit-learn** , with               RandomForestClassifier(n_estimators=2000,                            max_features="auto",                             max_depth=None)      After 70 minutes, I got a 2.9% test error, but with n_estimators=200 instead, I got a 2.8% test error after just 7 minutes. With **OpenCV** , I tried               rf.train(images.reshape(-1, 28**2),               cv2.CV_ROW_SAMPLE,               labels.astype('int'))      This ran for 6.5 minutes, and using `rf` for prediction gave a test error of 15%. I don't know how many trees it trained, as their Python binding for Random Forests seems to ignore the `params` argument, at least in version 2.3.1. I also couldn't figure out how to make it clear to OpenCV that I want to solve a classification problem, rather than regression -- I have my doubts, because replacing `astype('int')` with `astype('float32')` gives the same result. In **neural networks** , for the _permutation-invariant_ MNIST benchmark, the state of the art is 0.8% test error, although training would probably take more than 2 hours on one CPU. Is it possible to do much better than the 2.8% test error on MNIST using Random Forests? I thought that the general consensus was that Random Forests are usually at least as good as kernel SVMs, which I believe can get a 1.4% test error.


title: How can I generate a plot similar to that produced by plot.bugs and plot.jags from an mcmc.list?

text: R seems to be able to output nice summary plots from the `bugs` and `jags` objects generated by the functions R2WinBUGS::bugs and R2jags:jags. However, I am using the `rjags` package. When I try to plot the results of the function `rjags::coda.samples` using `R2WinBUGS::plot.mcmc.list` the results are diagnostic plots (parameter density, chain time series, autocorrelation) for each parameter. Below is the type of plot that I would like to produce, from Andrew Gelman's tutorial "Running WinBuugs and OpenBugs from R". These were produced by using the `plot.pugs`. The problem is that `plot.bugs` takes a `bugs` object as an argument, while `plot.mcmc.list` takes the output of `coda.samples`. Here is an example (from the `coda.samples`):                library(rjags)      data(LINE)      LINE$recompile()      LINE.out <- coda.samples(LINE, c("alpha","beta","sigma"), n.iter=1000)      plot(LINE.out)      What I need is either   * a way to generate a similar, information-rich, one-page summary plot similar to the one produced by `plot.bugs`   * a function that will convert `LINE.out` to a bugs object or ![enter image description here](http://i.stack.imgur.com/BSMmd.jpg)


title: Pointers to understand "rlm" in R better?

text: Could you please point me to tutorial/notes that can help me understand "rlm" better? Here is an example: `summary(rlm(stack.loss ~ ., stackloss, weights=myweights))` My questions are: a. How does the iterative procedure work?   b. What's the relation between the weights I supplied above "myweights" and the final used weights? My understanding is that the final weights used are actually the Huber loss weights? So are myweights not used at all? Ultimately, I want to use myweights... I have a feeling that actually I should put "myweights" into the "w" argument not the "weights" argument... c. Could you please give an example showing whether the final used weights are actually the Huber weights or my weights? I am also confused about the manual: The input arguments: > wt.method are the weights case weights (giving the relative importance of > case, so a weight of 2 means there are two of these) or the inverse of the > variances, so a weight of two means this error is half as variable? > > w (optional) initial down-weighting for each case. > > init (optional) initial values for the coefficients OR a method to find > initial values OR the result of a fit with a coef component. Known methods > are "ls" (the default) for an initial least-squares fit using weights > w*weights, and "lts" for an unweighted least-trimmed squares fit with 200 > samples. > > The returned values: > > w the weights used in the IWLS process > > wresid a working residual, weighted for "inv.var" weights only. Anybody please shed some light? Thank you!


title: Discrete distribution fitting

text: I have a discrete sample, numebr of passengers. when I try to find the best fit to this data set, none of the discrete distribution give the good fit according to goodness of fitness test. but the probability density function seems to be ok for example for Neg Binomial Should I rely on the plot or goodness of fitness test result? If I rely on GOF, all the distribution rejected. That means nothing fit to my data? I just want to use discrete distribution. and I have tried all the 8 famous discrete distributions. Here is my number of passanger data:                Passengers |      Freq.     Percent        Cum.     ------------+-----------------------------------               0 |         18        6.19        6.19               1 |         27        9.28       15.46               2 |         37       12.71       28.18               3 |         46       15.81       43.99               4 |         37       12.71       56.70               5 |         24        8.25       64.95               6 |         40       13.75       78.69               7 |         20        6.87       85.57               8 |         10        3.44       89.00               9 |         14        4.81       93.81              10 |          4        1.37       95.19              11 |          9        3.09       98.28              12 |          4        1.37       99.66              13 |          1        0.34      100.00     ------------+-----------------------------------           Total |        291      100.00


title: State space representation using KFAS package

text: I am using `KFAS` package for `R`. You can run               install.packages("KFAS")     library(KFAS)     ?regSSM      to see how this package allows to build a state space representation of linear regression models and many others. Now let we have the following state space system: $S_{t}=\alpha+(1+k_{t})L_{t}+v_{t}$ $k_{t}=\phi k_{t-1}+(1-\phi)\bar{k}+w_{t}$ being $\bar{k}$ a constant, a.k.a. the unconditional mean of the unobservable AR(1) process. Anyone can tell me how may I set this state space representation in `KFAS` through `regSSM` or any other `KFAS` package's function (like `arimaSSM`)?


title: How can I evaluate the accuracy of a clustering when I don't have information on the true class labels?

text: Already classified data set for the t-shirt factory problem I want to calculate the accuracy of my algorithm. I have the training data without any size information and I couldn't find the classified version of the data set. Should I create artificial data or find another data set?


title: Looking for catalog of "mechanism classes" that give rise to specific curve shapes

text: (I apologize for the length of this post. I don't know how to frame the question more succinctly.) I have some experimental data, in the form of a collection of curves with fairly little noise, which I'm trying to interpret. Specifically, I'm working on the question: what does these data say about the "internal workings" of the system that gave rise to them? After a bit of brainstorming, I came up with a couple of distinct candidate underlying processes/mechanisms, but I suspect that my brainstorming is only scratching the surface. I went online in search of either a more comprehensive collection of candidates, or else a general, abstract characterization of all such processes/mechanisms. Here's where I realized that I did not know how to search for what I was after. At this point, an example of what I'm after may help. Suppose that the shape of the curve of interest was a nice normal/Gaussian "bell". Among the many ways that one can characterize this curve is Galton's much-celebrated quincunx, which nicely demonstrates not only the normal bell shape, but its being the $N\to\infty$ limit of binomial distributions. So, if the data I were trying to interpret consisted of bell-shaped curves, the resource I'm searching for is one that would include the quincunx (or something close to it) as a prototype "underlying mechanism/process" for a normal bell-shaped curve. Now, given some other well-known curve shape (e.g. exponential decay, logistic curve, whatever), how does one search for candidate _underlying processes_ (analogous to the quincunx for the normal bell)? I'm having a really hard time searching for such information. I'm looking for any information that may help me with this sort of search. (The family of curves I'm interested in all lie in an envelope formed by a sigmoid-shaped curve with positive derivative everywhere and the constant curve corresponding to the asymptote of the sigmoid in the negative direction. The curves do not intersect each other.)


title: Reporting : Risk ratio vs Odd ratio

text: I have employee attrition model. I have reported odd ratios.My question is " can i use risk ratio?". Example : The odds of male to get attrited is 2.5 times than males. It is non-intutive for layman. It is better to say the chances of male to get attrited is more than males. I have googled and it yielded the relative risk is meaningless in case-control study. This term has origin in healthcare industry. I need to understand " Why researchers don't report relative risk ratio in HR,retail and telecom domains?". Should i report relative risk instead of odd ratio to make layman understand my findings? I know the difference between relative risk ratio and probability. I know how to calculate them. My question is more related to reporting of these statistics.


title: Differing p values (significant and insignificant) for the same predictor variable in different AIC models

text: I'm looking at a multinomial logistic regression analysis of deer behavioural responses to camera traps. The levels of the response variable are `no reaction`, `reaction` and `strong reaction`. I've selected a number of models based on their AIC values. However, I've found that the _p_ values gain and lose significance within differing models. For instance, with only `season` as a predictor:               Season               Summer  Winter          Coefficient r         0.8     0.03     Coefficient sr        0.4     0.5     P value     r         0.041   0.925     P value     sr        0.18*   0.024          AIC  1005.023     * P value of interest      With `species`, `camera model`, and `season` as predictors:               Species, Camera model and Season                              Muntjac  Roe   Sika  Bushnellb   Reconyx  Summer Winter          Coefficients  r       0.94   .43  -0.56    0.76        0.71    1.6    0.5     Coefficients  sr      1.11   .19  -0.19    0.28        0.77    0.95   0.8     P value       r       0.1    0.4   0.4     0.13        0.1     0.008  0.3     P value       sr      0.005  0.6   0.6     0.4         0.01    0.03*  0.02          AIC 1006.618     * P value of interest      Am I right to assume that this may be due to collinearity inflating the significance of this _p_ value, and therefore when analysing the probability of the predictor variable summer influencing the response variable `sr`, using the _p_ value from the first model? Is it correct to assume that the _p_ value from the 2nd model is due to collinearity and thus not representative of the actual significance of this variable (or level?) summer? Thanks!


title: Interpreting conflicting results from Random Forest & Logistic Regression?

text: I am using SKLearn and Statsmodel in python to build a RF and Logistic Regression, respectively. I have a feature that the RF indicates is important (feature importance of 0.202, closely behind #1 and #2 most important features). However, in running the logistic regression, the coefficient associated with this same feature is 0.0009, nearly 0. What is going on here? Why would splitting on this variable lead to higher information gain in the Random Forest, but contribute little to the logistic regression model?


title: Calculating sample size. Repeated measures within subject

text: I'm trying to calculate the sample size for a repeated measures within subject experiment. I'm going to test the subjects under 4 conditions and for each condition 6 measurements are taken at different time points. I want to use an alpha level of 0.05, but I'm unsure how to calculate the sample size needed. I want to find the sample size required to see a significant difference, between each condition, if the 6 time points are averaged. I also want to find the sample size required to see a significant difference, when comparing for example condition 1, time point 2 with condition 2, time point 2 and so on. Can anyone help me? I'm using G*Power 3. So far I've tried using the following ![G*power 3, values](http://i.stack.imgur.com/SMGkZ.png)


title: Goodness of fit (Kolmogorov-Smirnov) test for a time series?

text: This is a somewhat basic question, I guess. I have a sequence of random variables $X_1,\ldots,X_n$ that I believe to be i.i.d. uniform on $[0,1]$. Being i.i.d. uniform corresponds to my model predicting the correct distribution for $P(Y_k\mid Y_{1:k-1})$ for some time series $Y$. I can test that by, say, applying the Kolmogorov-Smirnov or the Anderson-Darling test. However, there is a type of dependence that I want to check: all together they might have a uniform ecdf, but, $X_1,\ldots,X_{n/2}$ and $X_{n/2+1},\ldots,X_n$ separately might not. But if I apply the KS test to each half, the result will not be independent from the statistic generated by the KS test for the whole sequence. Furthermore, if I apply the KS test to each window $X_k,\ldots,X_{k+m}$, and get a function $k\mapsto \text{p-value}$, $k=1,\ldots,n$, I can check visually whether there are intervals where KS says they are not uniform. But there are $n-m$ tests for each $k$, so they are very much not independent of each other, so this isn't very sound. What goodness-of-fit test should I use here? How do I justify this procedure properly? Should I be doing this differently? Which textbooks/articles do I need to read?


title: R: Setting an F-statistic to determine variables for a multiple linear regression model

text: I am trying to understand the steps behind the linear regression process. I already have a linear model like: `lmodel1 <\- lm(y~x1+x2+x3, data=dataset)` for which R calculates several different things (`Coefficients, Intercept, Residuals, F-statistic` and `p-value`) among others. At this point, i am mostly interested in `F-statistic` and `p-value`. So far, i have concluded to the following: The process is iterative and begins taking under consideration every variable. In order to achieve an optimal `y` some `x` variables have to be "taken out". This comes as a result of calculating F-statistic, which quantifies the importance of each `xi` and the dependent variable `y`. When `F value` is smaller than `p-value`(?) that variable is removed. Next step of the process is to compare that `F-statistic` of a `xi` independent variable, with an `F-to-enter` and `F-to-remove` in order see if the removed variable will be re-inserted to the equation.(?) Now, please do correct me if i am wrong regarding the steps desribed above. Is that what happens under `lm()`'s hood. Are those the right variables.? R-wise speaking how does these values can be shown, inserted or calculated in a multiple linear regression model.? How is `ANOVA` related to the above? I am afraid R's `summary` and `help` take too much for granted. Thanks in advance for any suggestions.


title: Expression for conditional density for ARCH processes

text: I am reading Stephen Taylor's Asset Dynamics book and came across something I didn't fully understand. For an ARCH process, the return series is modeled as $r_t = \mu_t + h_t^{1/2}z_t$ where is $z_t$ is $D(0,1)$ and may be normal/non-normal and $h_t$ is conditional variance (some function of subset of $\theta$). Then he says the conditional density of $r_t$ is $f(r_t|\theta) = \frac{f(z_t(\theta)|\theta)}{\sqrt{h_t(\theta)}}$. I don't understand how does the conditional variance term occur in the denominator of density relation and why doesn't $\mu_t$ occur in the relation. Can someone help me how to derive this?


title: Can you express a probability density as a percentage?

text: I'm working on a naive Bayes classifier that calculates probabilities using a normal Gaussian distribution. This works very well when I am classifying something into two mutually exclusive buckets (e.g. spam vs. not-spam), but when I am working with a factor that is not easily classified that way (when the classifications are not mutually exclusive) I would like to express the result as a percentage. When I combine the probability density of several factors (by multiplying them together) I tend to get a very small number and I would like to adjust that so I can express it in a 0-100 percent range, so it will be more easily understood. Is there another factor I can use to adjust the probability density into a percentage range? For example: in the Wikipedia article for naive Bayes classifiers, there is an example of using height, weight and shoe size variables to classify a person as male or female. After computing the numbers, the posterior result for the female classification is 5.3778E-4 (or .00053778). Out of context that seems minuscule, but if you compare that number to the result for the male classification, the percentage would be 99.99% female. What factors (if any) could I apply to the posterior result to get that percentage, without the benefit of the male result to compare it to?


title: What is the difference between Z-scores and p-values?

text: In network motif algorithms, it seems quite common to return both a p-value and a Z-score for a statistic: "Input network contains X copies of subgraph G". A subgraph is deemed a motif if it satisfies   * p-value < A,   * Z-score > B and   * X > C, for some user-defined (or community-defined) A, B and C. This motivates the question: > **Question** : What are the differences between p-value and Z-score? And the subquestion: > **Question** : Are there situations where the p-value and Z-score of the > same statistic might suggest opposite hypotheses? Are the first and second > conditions listed above essentially the same?


title: Testing for a drop in bookings

text: We're developing real-time alerts for fine-grained (every 5 minutes) time series bookings data, and I'm looking for the best approach to doing this. Idea is that if over the past 10–15 minutes (say) there's a big drop in bookings volume relative to some expectation, then I want to fire an alert. We've been using historical data to set the expectation. I'm familiar with time series analysis and exponential smoothing. Holt–Winters came up as a possibility here, but even triple exponential smoothing appears to handle only one seasonality, and bookings data in my industry (travel) definitely has at least three driving seasonalities: daily, weekly and annual. **What we've been doing.** One approach we've been using – and it seems to work OK in most cases – is to do a difference between two sample means hypothesis test. Say that we break bookings up into five-minute blocks. Then we build a size three sample over the past 15 minutes, a size 12 sample for the corresponding 15 minute spans over the past four weeks, and then test the difference in their means. I am a concerned that the sample sizes (especially for the current bookings sample) are too small to be trustworthy. The test accounts for sample size, and anyway the variance within a sample is generally quite small, at least at volume, so maybe this is OK. But using an $n=3$ sample feels wrong. Besides this approach neglects the annual cycle, which I'd like to be able to capture if possible. **Extra background.** A couple of considerations to add. Not sure about their relevance, but just in case: 1) These are global bookings (i.e. people all over the world book), so the bookings volumes are not simple sinusoids. They are roughly the sum of multiple bookings sinusoids corresponding to different regions (NA/LATAM, EU, APAC) in addition to those associated with the different seasonalities. I know we can analyze the individual regions to tease apart the sinusoids, and we want to do that. But also we want to know when the aggregate volume has dropped. 2) There is some sense in which bookings seem to have a binomial distribution. If we fix the user volume, then we can think of a booking as a "success" in a binomial trial. This comes up for example if you have a web server farm with many servers, and you look at the distribution of bookings counts across the servers within a short time window (e.g. 30 minutes). That distribution passes a normality test (I've been using Shapiro), and I assume it's because it's a binomial distribution given the underlying mechanism. Bookings over time aren't binomially distributed because the user volumes themselves fluctuate widely. But I still wonder if somehow we can use the large number of underlying binomial "conversion opportunities" as a way to drive an appropriate test. My question is: what approach should we be using to test for bookings drops, given the above? **Update.** Based on what Stephan says below (that this is really anomaly detection, not forecasting), I have another idea. What if I were to build a comparison sample based on historical 5-minute bookings counts (e.g. seven 5-minute blocks from the previous week, flanking and including the block corresponding to the current one, another seven from the week before that, and so on, maybe four weeks back, giving me an $n=28$ sample), and then use an interquartile range-based outlier test? Simple, and the detection bands are given by the boxplot whiskers. Is this at all plausible?


title: How to show that an estimator is consistent?

text: Is it to show that MSE = 0 as $n\rightarrow\infty$? I also read in my notes something about plim. How do I find plim and use it to show that the estimator is consistent?


title: Principles to create natural data for a clustering algorithm?

text: I have come up with an algorithm to cluster geospatial data points. I have quite a few volunteers (100) collecting data for me, using my app on their smart phones to check in at places. However, I'm afraid the data will be positively biased since the volunteers all receive instructions on how to use the app. The data will just be too perfect and easy for my algorithm. So, I intend to ask some to deliberately create noisy data (by checking in at wrong positions). My questions are:   * Is asking people to create wrong data a sound way?    * If so, what is the optimal percentage of such data? Say, 20 out of 100?  Better yet, can you point me to some good research papers on **principles** of creating effective data? I.e, the data that best replicates real-world situations. Thanks.


title: Logistic mixed model

text: In the logistic mixed model ${\rm logit}(P(Y_i=1))= α + βX_i + u_i + ε_i , i=1,...,m$, when we know $u_i\sim \mathcal N(0,σu^2)$, and $ε_i\sim\mathcal N(0,σi^2)$, and if we know $σi^2$ in each area $i$:   1. How do we use $σi^2$ in estimation of parameters (likelihood formula)?   2. How can we estimate the $σu^2$?


title: What is a good technique for grouping objects based on binary or dichotomous traits?

text: I have a set of objects each of which has a list of traits. Data on the traits is binary: an object has a trait or does not. The number of objects that I have is moderately greater than the number of traits, and much greater if you exclude traits that are held by only a few objects. Virtually all objects have at least 3 or four traits, and often as many as 15. Currently my sample has about 200 objects in it, but I may find more I am looking for sensible candidate algorithms for dividing the objects into groups based on these trait measurements. In the best of all possible worlds, the objects would divide neatly into distinct groups that share traits with other group members but not with outsiders. I already know my data is not so neat. I am looking for an algorithm that assigns objects to groups in such a way that the number of groups is determined by the data, or the data plus a small number of cutoff parameters, and for which a grouping is preferred based on some combination of the following three properties:   1. the sharing of traits within a group is **high** ;    2. the sharing of traits between members of different groups is **low** ; and    3. for the members of each group, the pattern of sharing traits between groups is **similar** , e.g. if member 1 of group A shares 50 percent of its traits with members of group B but only 1 percent with members of group F, then members 2 through 20 of A will also have high rates of sharing with group B and low rates with group F. Is there a standard or a short list of common algorithms that have been used to address problems with this structure? Alternatively, is there an algorithm that someone would like to champion for this problem? I am hoping to find an algorithm that I can implement with R (or with a spreadsheet). I do not know the relative importance of criteria 1, 2, and 3 and would ideally like to be able to vary them.


title: Standard Deviation Grows Quadratically with Input Variable

text: I take an input `x`, based on which I do an experiment that gives me several data points. I compute the standard deviation of these data points. Then, I change `x`, and repeat the process. Finally, I ask the dependence of the standard deviation on `x`. I am observing a quadratic dependence on `x`. Is there anything that can be inferred about the process that is picking the data points? Perhaps it may be helpful to briefly mention the origin of this problem – I'm measuring noise in an experiment and trying to characterize the origin of the noise by characterizing its effects. `x` is a variable experimental setup, and the process that picks my data points is the noise. As a function of `x`, the noise presumably picks the data points in such a way that the distribution depends quadratically on `x`.


title: Relationship between regressing Y on X, and X on Y in logistic regression

text: Correlation and linear regression are sometimes distinguished in statistics books by saying that the former is symmetric and the latter is asymmetric in the following sense: in the case of correlation, no distinction is made between dependent and independent variables, whereas it makes a difference which variables are treated as dependent and independent variables in a regression equation. To be sure, it is possible to treat $Y$ as a predictor of $X$ instead of treating $X$ as a predictor of $Y$ by using the following equations: $X$ as a predictor of $Y$:   \begin{align} \text{Slope: }b_1 &= r_{XY} \cdot \frac{(sY)}{(sX)} \\ ~\\ ~\\ \text{Intercept: }b_0 &= \bar{y} - b_1 \bar{x} \end{align} $Y$ as a predictor of $X$:   \begin{align} \text{Slope: }b_1 &= r_{XY} \cdot \frac{(sX)}{(sY)} \\ ~\\ ~\\ \text{Intercept: }b_0 &= \bar{x} - b_1 \bar{y} \end{align} However, the regression lines to which the scatter plot will be fitted will differ depending on whether $X$ is treated as a predictor of $Y$ or $Y$ is treated as a predictor of $X$. My question is whether something similar holds in the case of logistic regression and whether it is possible to formulate similar equations for the regression coefficients in logistic regression for when $Y$ is used as a predictor of $X$.


title: Is there a better way than side-by-side barplots to compare binned data from different series

text: I have a number histogram from a data source. I want to compare their distribution. For example series 2 is more to the right than series 1. ![enter image description here](http://i.stack.imgur.com/7yF5L.png) Is there any good way to contrast besides plotting the histogram side by side. I can think of qqplot or box plot. But the data is already binned. Would it make sense if I interpolate the percentile linearly? EDIT1: I am exploring the chart in Excel. There not a lot of flexibility for creativity. I hacked Excel to make it in demographic chart style. The good thing is it is a familiar visualization to many people. ![enter image description here](http://i.stack.imgur.com/la4dw.png) EDIT2: I interpolate the data to get 25 and 75 percentile to make this interquartile, not quite boxplot, chart. It loses a lot of details. But it may be useful for comparing a lot of series. ![enter image description here](http://i.stack.imgur.com/au5KC.png) EDIT3: Sorry for not posting the raw data initially. The stuff is disorganized while I'm experimenting.               16-19,0.0028826,0.031066407     20-24,0.05581761,0.094111701     25-34,0.278301887,0.215492493     35-44,0.289046122,0.214615109     45-54,0.228773585,0.230744934     55-64,0.120807128,0.163116187     65+,0.024633124,0.050853168


title: What Machine Learning Algorithm would be appropriate for this scenario

text: I have a PHP/MySQL application that stores symptoms and the appropriate drug. What machine learning algorithm should I use to predict the drug for any symptoms. Also, what would be the format of the training set? Features: 2 (not numeric, string/character) Samples: More than 100 No Missing Values The Labelled data would be as follows: A combination of symptoms relate to a combination of drugs. So, (symptom1, symptom2, symptom3) --> (drug1, drug2) (symptom1, symptom2, symptom4) --> (drug2, drug3, drug4) In general, one ( or more than one ) symptom links to one ( or more than one ) drug: (symptom1, etc...) --> (drug1, etc ...)


title: Control chart Xbar and R missing values

text: I am taking a business statistics course and I can't seem to figure out how to find the missing figures in the following problem; here's the question: Using the data below, determine the trial central line and trial control limits for the Xbar chart. Subgroup size is 5. For the R chart, the trial central line is 0.11, the UCL is 0.23, and LCL is 0. Note not all of the subgroups are included in the data; however, the sums are correct. (5 subgroups)                   X-Bar    R     1   2.08    .05     2   2.09    .24     3   2.08    .10     *   *       *     *   *       *     *   *       *     25  2.09    .11         =52.0   =2.75


title: Using EM algorithm for record linking

text: I am interested in linking records across 2 datasets by first name, last name, and birth year. Might this be doable with the EM algorithm, and if so, how? Consider the following record in the 1st as an example: Carl McCarthy,1967. I will search through all records in the 2nd dataset, and assign a jaro-winkler distance between the 1st name and Carl and a jaro-winkler distance between the last name and McCarthy. These distance are probabilistic as is the distance between the birth years. We combine those 3 probabilities (multiply? average?) into 1. Now comes the decision rule part. Let us rank all of the probabilities from highest to lowest. First, we want P(first hit is match) >= threshold. Second, we also want P(first hit is match) / P(second hit is match) >= threshold if P(second hit is match) exists. Third, we want the first hit in this second dataset to match for no more than 1 person in the 1st dataset with Carl McCarthy,1967. How may these thresholds be determined? I prefer approaches in Stata and/or Perl. See, for example: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1479910/pdf/amia2003_0259.pdf (Although with that, I still do not fully follow the why or how, and what the inputs and outputs are, as well as the assumptions and how restrictive they are).


title: Modeling threshold effects

text: First, let me give some context for my question. I'm working on a project in which our research team implemented a workshop intervention at a number of different departments at my university. Each department that received the intervention also had a matched control department. We obtained several different measures from professors and staff in these departments at three timepoints: before the intervention and at two follow-up time points. To examine the effects of our intervention on our various outcome measures, we are using the following model using the lme4 package, where pair is a code indicating the department pair the respondent was from, dept is a code indicating the department the respondent was in, and person is a code for each respondent: mod <- lmer(outcome ~ condition * time + (1|pair) + (1|dept) + (1|person), data = d) We also have records on the percent of people that attended the workshop within the intervention departments. Using these attendance data, we could fit the following model: mod <- lmer(outcome ~ attendance * time + (1|pair) + (1|dept) + (1|person), data = d) However, we would like to go one step further by estimating the attendance threshold required for the intervention to affect the outcome. I've reached a bit of a blind alley in terms of how to test for threshold effects. Does anyone have any ideas on how we might go about this?


title: Feature selection with Caret for data with more than one target

text: I am trying to do some feature selection, having around 3500 variables for about 200 samples. To each sample is associated two numerical values (the expected outcome). I can't manage to make the caret work with this, or even find any information on this. Does anybody know how to do this? As an example, my data is roughly in the following format: Samples:               S1    2.1 1.2 3.1 ... 4.2 1.7 5.2     S2    3.4 1.1 4.5 ... 5.3 1.2 5.7     ...     S3499 2.4 3.5 5.1 ... 2.2 1.5 5.7     S3500 4.1 1.2 5.4 ... 1.2 2.1 5.8      Targets:               S1    1.82 1.44     S2    2.44 1.22     ...     S3499 1.23 1.32     S3500 1.99 1.51      Thanks, Swatchpuppy


title: ARIMA forecasting

text: I have a data that I would like to use arima model to perform forecast. when I use auto.arima, my results does not seem right. When I change my arima order to c(1,0,1), numbers starts changing and seem more realistics. How do I know which arima order to pick? Which orders set I should try to compare results? Not sure but sigma seems like the error. Maybe I should pick the order level that provides the smallest sigma? Any guidance is greatly appreciated. Here is my data set: tmp                         Jan      Feb      Mar      Apr      May      Jun      Jul      Aug      Sep      Oct      Nov      Dec     2011                                              38.96315 38.27304 40.79011 40.59401 42.22266 42.76763 38.69254     2012 42.77221 45.47796 44.32276 46.30973 51.64651 47.91630 49.75013 24.68711 24.39294 24.35242 27.60628        ar<-auto.arima(tmp) ar result:               Series: tmp      ARIMA(0,1,0)                              sigma^2 estimated as 43.54:  log likelihood=-56.2     AIC=114.4   AICc=114.66   BIC=115.23      perform forecast on ar:               vol_fore <- forecast(ar,h=12)      result of vol_fore:                       Point Forecast      Lo 80    Hi 80      Lo 95    Hi 95     Dec 2012       27.60628 19.1496201 36.06294  14.672936 40.53962     Jan 2013       27.60628 15.6467578 39.56580   9.315770 45.89679     Feb 2013       27.60628 12.9589169 42.25364   5.205072 50.00748     Mar 2013       27.60628 10.6929623 44.51959   1.739594 53.47296     Apr 2013       27.60628  8.6966162 46.51594  -1.313554 56.52611     May 2013       27.60628  6.8917813 48.32077  -4.073811 59.28637     Jun 2013       27.60628  5.2320644 49.98049  -6.612129 61.82468     Jul 2013       27.60628  3.6872375 51.52532  -8.974738 64.18729     Aug 2013       27.60628  2.2363044 52.97625 -11.193748 66.40630     Sep 2013       27.60628  0.8639778 54.34858 -13.292541 68.50510     Oct 2013       27.60628 -0.4412831 55.65384 -15.288765 70.50132     Nov 2013       27.60628 -1.6884441 56.90100 -17.196133 72.40869      forecasted data points are the same. Therefore, I dont think auto.arima is a good model for this. I've been looking at order levels, but could not figure out which ones to pick.


title: Bias in Estimators of Lognormal

text: I am modelling a process distributed as a 2 parameter lognormal distribution; determining the parameters by maximum likelihood. I have simulated the bias in the estimators (logmean and logsd) as well as the skewness. The logmean seems unbiased and the logsd and skewness seem asymptotically biased; underestimating these parameters. Is there any closed form correction of the bias of these estimators?


title: Variable reduction techniques on qualitative factor with several variables

text: Suppose that I want to regress the income level of a worker with qualitative variables such as eye color and the state that they live. Obviously, if I try to create a model with these variables, I will have multiple categorical variables, such as a variable representing a person living in California or having green eyes. Is there a way to reduce the amount of categorical variables and combining these categorical variables into a new category? For instance, I want to find a way to categorize the different states into 3 groups and I do not want the eye category to be in the group. Are there techniques such as PCA and factor analysis that can do that? Would you use discriminant analysis for that technique? I know I can categorize the United States by North, South, East and West, but I do not want to do that. I feel that is a little too lame. Please let me know if I should clarify anything.


title: Study that uses efficient sample size

text: Suppose we have $10$ pairs of matched data (i.e. pre-exposure and post exposure). Call the exposure $E_1$. So the data are in the form $(\overline{x}_1, \overline{y}_1), \dots, (\overline{x}_n,\overline{y}_n)$. For example, $(\overline{x}_1, \overline{y}_1)$ is the average of measuring $E_1$ of subject $1$ 10 times pre-exposure and $10$ times post-exposure. Suppose we have some other exposure $E_2$. We want to test whether $E_2$ is effective by itself. We also want to compare $E_2$ and $E_1$. We only have the data on $E_1$. Also we are primarily interested in $E_2$. What are some ways of reducing the sample size required for these experiments? A crossover study seems good because each person can serve as their own control. Also the sample size for both problems would be the same.


title: Factor analysis using stata "predict" command and get negative value for non-negative variable?

text: I am running factor analysis on stata to reduce a few variables to a single explanatory variable which means "experience" of a manager (should be non- negative value), however, after using "predict" command I check the range of the new variable and found that there are many negative values, how do I avoid this? The original variables are 1) number of projects done by a manager 2) number of people interacted for each project 3) number of unfinished projects. Since these 3 variables are highly correlated so I want to use factor analysis to reduce them to one single variable. If there is any other method that can also tackle this problem, please let me know! Thank you guys!!


title: MLR prediction intervals joint region

text: This thread offers some good information on the calculation of linear regression prediction intervals and includes a link to some practical notes. At the end of the link the author claims that, in case of multiple linear regression, we should consider the joint region of the predictor variables. How should MLR prediction intervals be interpreted in light of their joint region limitations?


title: How to split r-squared between predictor variables in multiple regression?

text: I have just read a paper in which the authors carried out a multiple regression with two predictors. The overall r-squared value was 0.65. They provided a table which split the r-squared between the two predictors. The table looked like this:                           rsquared beta    df pvalue     whole model     0.65   NA  2, 9  0.008     predictor 1     0.38 1.01 1, 10  0.002     predictor 2     0.27 0.65 1, 10  0.030      In this model, ran in `R` using the `mtcars` dataset, the overall r-squared value is 0.76.               summary(lm(mpg ~ drat + wt, mtcars))          Call:     lm(formula = mpg ~ drat + wt, data = mtcars)          Residuals:         Min      1Q  Median      3Q     Max      -5.4159 -2.0452  0.0136  1.7704  6.7466           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   30.290      7.318   4.139 0.000274 ***     drat           1.442      1.459   0.989 0.330854         wt            -4.783      0.797  -6.001 1.59e-06 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 3.047 on 29 degrees of freedom     Multiple R-squared:  0.7609,    Adjusted R-squared:  0.7444      F-statistic: 46.14 on 2 and 29 DF,  p-value: 9.761e-10      How can I split the r-squared value between the two predictor variables?


title: Plot log-normal distribution in R

text: I need to plot lognormal distribution with mean 1 and variance 0.6 in `R`. I tried to do this using `rlnorm` function in `R` as               x= rlnorm(500, log(1), log(0.6))          plot(density(x))      `log(0.6)` is negative that might be the reason my code is not working but the `R` documentation for `lnorm` says the value for standard deviation is given in log. Correct me if I have understood it incorrectly. Can someone please give me correct `R` code?


title: Validity of grand mean intra-rater Kappa or ICC

text: For 5 raters, I have calculated intra-rater Cohen's Kappa statistics for the test-retest nominal ratings, intra-rater Kendall's tau-b statistics for the test-retest ordinal ratings, and intra-rater ICCs for the test-retest scale ratings. Is it valid to calculate overall mean intra-rater reliability coefficients (Kappa, Kendall's tau-b, ICC, as appropriate) across all raters for each of these variables? If so, can I simply average the individual raters' reliability coefficients for overall means?


title: Does finite kth moment imply lesser moments are finite?

text: > **Possible Duplicate:**   >  Proof that if higher moment exists then lower moment also exists For a random variable $X$, lets say I know $E[X^k]$ is finite and I know that $E[X]$ is finite. Can I say that all moments between first and kth moment are finite?


title: Sample variance order

text: Is it true that (and if so, how does one prove) the following. $$E\left|\hat{Var}_{n}(X)-Var(X)\right|^{2}=O(n^{-1})$$ where: • $X$ is a random variable with mean $\mu$ and variance $\sigma^{2}$ • $\hat{Var}_{n}(X)$ is the sample variance of $X$ from $n$ i.i.d. random variables $X_{1},\cdots,X_{n}$ with mean $\mu$ and variance $\sigma^{2}$. Many thanks in advance. (Feel free to change my notation).


title: Why does a 95% CI not imply a 95% chance of containing the mean?

text: It seems that through various related questions here, there is consensus that the "95%" part of what we call a "95% confidence interval" refers to the fact that if we were to exactly replicate our sampling and CI-computation procedures many times, 95% of thusly computed CIs would contain the population mean. It also seems to be the consensus that this definition does _not_ permit one to conclude from a single 95%CI that there is a 95% chance that the mean falls somewhere within the CI. However, I don't understand how the former doesn't imply the latter insofar as, having imagined many CIs 95% of which contain the population mean, shouldn't our uncertainty (with regards to whether our actually-computed CI contains the population mean or not) force us to use the base-rate of the imagined cases (95%) as our estimate of the probability that our actual case contains the CI? I've seen posts argue along the lines of "the actually-computed CI either contains the population mean or it doesn't, so its probability is either 1 or 0", but this seems to imply a strange definition of probability that is dependent on unknown states (i.e. a friend flips fair coin, hides the result, and I am disallowed from saying there is a 50% chance that it's heads). Surely I'm wrong, but I don't see where my logic has gone awry...


title: What statistic method using spss to use to analyze online buying behavior based on different variables?

text: I have the data with me however I don't have any idea as to how to analyze the variables. Can anybody help me? The variables/questions are as listed:               1. Gender        2. Age     3. Education Level     4. Occupation     5. Household monthly income (in PHP $)       6. Daily usage of the Internet       7. Have you ever used the internet for purchase? (If your answer is YES, please skip questions 14-20. If your answer is NO, please skip questions 8-13 and jump to the next page)        8. How many times have you purchased online?              "Reasons for adopting the online purchases          9. Lower prices     10. Easiness of online buying procedures     11. Wide variety of products         12. Various payment options      13. High quality of products         14. Security and privacy reasons         15. Need to physically examine the product       16. Prefer to buy from retail stores         17. Do not own a credit card              "Reasons for not adopting the online purchases     18. Are unaware of the buying procedure through the Internet         19. Shipping delays          20. Unaffordable shipping fees       General expectations:     21. Online stores should have a good reputation      22. Online stores should provide adequate payment options.       23. Delivery service should be reliable          24. The internet consumers must be protected by government laws.         25. Consumers can return a product and get refunded the purchase price       26. Security mechanisms are used to ensure personal data safety      27. The company has also retail stores       28. Online stores offer guarantee for their provided products        29. Shop any time of the day         30. Shop abroad      31. Saves time       32. Easy to do comparison shopping between products, as well as, online stores       33. More easily you can find a product compared to retail stores         34. Have more options compared to retail stores          35. Easy to find real bargains       36. Provided products are cheaper compared to retail stores      37. The whole buying procedure is easier compared to retail stores       38. Consumers can find products that are not in retail stores        39. Have much more time to evaluate and select a product         40. Online stores promise more than they can practically offer       41. Consumers can not completely trust online stores         42. Online stores are not always official representatives of their offered products      43. Consumers find it difficult to confirm the reliability of the provided products      44. It is possible to have your credit card data intercepted         45. It is difficult to have after-sales services          frequency of use for thee purposes:      46. Surfing (e.g. read news/articles)        47. Social Networking        48. E-mail       49. Job search       50. Research         51. Chatting (e.g. skype/yahoo messenger)        52. Playing Online Games         53. Searching for product info


title: Hessian of logistic function

text: I have difficulty to derive the Hessian of the objective function, $l(\theta)$, in logistic regression where $l(\theta)$ is: $$ l(\theta)=\sum_{i=1}^{m} \left[y_{i} \log(h_\theta(x_{i})) + (1- y_{i}) \log (1 - h_\theta(x_{i}))\right] $$ $h_\theta(x)$ is a logistic function. The Hessian is $X^T D X$. I tried to derive it by calculating $\frac{\partial^2 l(\theta)}{\partial \theta_i \partial \theta_j}$, but then it wasn't obvious to me how to get to the matrix notation from $\frac{\partial^2 l(\theta)}{\partial \theta_i \partial \theta_j}$. Does anybody know any clean and easy way of deriving $X^T D X$?