title: How to handle online time series forecast?

text: I have been dealing with the following problem. I have kind of a real time system and every time frame I read its current value, creating a time series (such as 1, 12, 2, 3, 5, 9, 1, ...). I'd like to know methods (statistical and machine learning) to forecast the next value in an online fashion (meaning every time a new value is read). I tried R's ARIMA and Weka's SMOreg, they result in good predictions, but they are kind of slow as every time a new value comes I have to redo the math. P.S. It would be great also if the method had a confidence interval.


title: Finding mean and standard deviation

text: Experience raising Labrador puppies revealed that the weight of the puppies at six months follow a normal distribution. The subsequent weights (in pounds) of a small sample of seven six-month puppies were: $$ \begin{array}{|c|c|} \hline \mbox{Puppy #}& \mbox{Weight} \\ \hline 1 & 16 \\ 2 & 46 \\ 3 & 34 \\ 4 & 31 \\ 5 & 21 \\ 6 & 26 \\ 7 & 22 \\ \hline \mbox{Sum} & 196 \\ \hline \end{array} $$ Based on this information: > a) What is the mean weight of these seven puppies? Which I got to be $28$. > b) Based on this small sample, what is the standard deviation? I got this to be $9.27$. > c) Construct the 95% confidence interval for the mean weight of all six- > month puppies based on this information. I got this to be $(9.83, 46.17)$ but I was wondering could I write this as $9.83 < \mu < 46.17$? > d) Would it be reasonable to state that this litter contained two “unusual” > puppies; i.e., a runt and a giant? Why or why not? I have no clue what this is asking. What I was thinking was we could say no because all of our data falls in our CI. Is this a good way to approach it?


title: Interpreting interaction terms with a log transformed dependent variable

text: In my analysis I'm using a linear regression with a log transformed dependent variable. I've included an interaction term in my regression model as I'm testing for interaction between particular independent variables. The model is as follows: ln(Y) = α + βX1 + βX2 + βX3 + βX4 + βX5 + β(X1*X2) with coefficients ln(Y) = 6.703 + -1.164X1 + 0.140X2 + 0.365X3 + -0.421X4 + 0.071X5 + 0.687(X1*X2) I understand that when you interpret the coefficients of the independent variables in the model, the coefficients refer to multiplicative changes in the dependent variable (i.e. X1's coefficient of -1.164 refers to [e(-1.164)-1]*100 = a decrease of 68% in the dependent variable when X1 is 1 given X1 is a dummy variable). What I don't understand is how to interpret the coefficient of the interaction term when there is a log transformed dependent variable. Using the coefficients above, should it be a 15% increase ([e(0.140)-1]*100) in the dependent variable for a 1 unit increase in X2 when X1 is 0 and a 128% increase ([e(0.140+0.687)-1]*100) in the dependent variable for a 1 unit increase in X2 when X1 is 1? Any guidance would be great!


title: Proper scoring rules for observations with different supports

text: Suppose to have a bivariate variable $z_t=(x_t, y_t)$ indexed by $t=1,2, ..., T$. Suppose now that the two components have different support, i.e. in my specific problem $x_t \in \mathcal{S}$, where $\mathcal{S}$ is a circle, and $y_t \in \mathcal{R}$. I estimate a model (Bayesian hidden Markov model) and then i want to evaluate how the model fit the data. In other situation, where i had only a univariate variable, i used the scoring rules: if my univariate variable was circular i used the continuous ranked probability score for circular variable, while when the univariate variable was linear i used the predictive log scoring loss. I am quite sure i cannot use a single scoring rules for bivariate distribution with different support. Am I right? What i need to find is a number that summarize the quality of forecast. Any suggestion?


title: Is ANOVA being misused for these binary data?

text: Everywhere I have read that if you want to test the hypothesis of dependence of a binary outcome on your data, you should not apply ANOVA (if its incorrect, please let me know). We can apply logistic regression, chi square testing etc. I have come across a paper in the reputed **Journal of Vision** , where I believe the ANOVA is being misused. I would like to mention that this is not just one case and I have seen many papers like this: Paper In this paper the authors are mainly seeing how inversion of face images leads to accuracy changes. However, they have written this statement: > ANOVA was conducted with the main factors of Orientation (upright or > inverted), Change Type (configural or featural) and Region (eyes or mouth). > The results (Figure 3) showed a significant main effect of Orientation, F(1, > 21) = 24.35, p < 0.01, and Region, F(1, 21) = 16.97, p < 0.01.` Can someone help me out by figuring how ANOVA is being applied here? All the 3 parameters Orientation, Change Type and Region are binary. However the authors are reporting the p values and hypothesis significance. Is this the correct way to do so? I am mainly asking this not just from the viewpoint of this paper. As I have mentioned many authors just take factors that are binary, apply ANOVA and show that one of them is better in a specific task.


title: Analysing a grouped 0 to 10 scale using ordinal logistic regression

text: I read in an article that the logit link is considered suitable for analyzing ordered categorical data evenly distributed among all categories. I want to do ordinal logistic regression and I have an ordinal variable on the scale of 0-10 and have to make three groups, 0-6, 7-8 and 9-10. Does this mean that the data is not evenly distributed and I shouldn't use ordinal regression or do they mean that there has to be approximately the same number of observations in each category?


title: When is it appropriate to select models by minimising the AIC?

text: It is well established, at least among statisticians of some higher calibre, that models with the values of the AIC statistic within a certain threshold of the minimum value should be considered as appropriate as the model minimizing the AIC statistic. For example, in [1, p.221] we find > Then models with small GCV or AIC would be considered best. Of course one > should not just blindly minimize GCV or AIC. Rather, all models with > reasonably small GCV or AIC values should be considered as potentially > appropriate and evaluated according to their simplicity and scientific > relevance. Similarly, in [2, p.144] we have > It has been suggested (Duong, 1984) that models with AIC values within c of > the minimum value should be considered competitive (with c=2 as a typical > value). Selection from among the competitive models can then be based on > such factors as whiteness of the residuals (Section 5.3) and model > simplicity. References:   1. Ruppert, D.; Wand, M. P. & Carrol, R. J. Semiparametric Regression, Cambridge University Press, 2003    2. Brockwell, P. J. & Davis, R. A. Introduction to time-series and forecasting, John Wiley & Sons, 1996 So given the above, which of the two models below should be preferred?               print( lh300 <- arima(lh, order=c(3,0,0)) )     print( lh100 <- arima(lh, order=c(1,0,0)) )      More generally, when is it appropriate to select models by blindly minimizing the AIC or related statistic?


title: Parameter Estimation

text: I have the data in the form of $Y \in \mathbb{R}^T$ a time series. For each point in time I have $ m $ real features $ f_i \in \mathbb{R}^m$. I want to use the following model to fit the data $ \hat{y}_t = \langle (y_{t-1}, ..., y_{t-h}), W_u(f_t)\rangle + B_v(f_t)$ where $ W $ stands for weight and $ B $ stands for bias. A parameterized linear model Both $W$ and $B$ are differential and chosen from a family of predefined functions (i.e. linear, polynomial, ...) mapping features to weight of previous values and bias. The parameterization of $ W $ and $ B $ which is $u, v$ is the desired goal. The loss function is the MSE of $ y_t - \hat{y}_t $. What methods could help me find $u$ and $v$ ? examples: $W:R^m \rightarrow R^h$ such that $ W_i(f) = \sum\limits_{j=1}^m u_{i,j}f_j$ or $ W_i(f) = \sum\limits_{j=1}^m e^{\frac{1}{u_{i,j}}}f_j$ or $ W_i(f) = \sum\limits_{j=1}^m {f_j}^{u_{i,j}}$ or $ W_i(f) = \sum\limits_{j=1}^m cos({u_{i,j}{f_j})}$


title: Generating and analyzing predicted data in JAGS

text: I am trying to model a logistic regression on data predicted from the results of a previous logistic regression, and am having trouble figuring out how to do it, either in OpenBUGS (BRugs) or JAGS (with rjags). I'm attempting to determine the probabilistic maturation reaction norms (PMRNs) of different fish populations. Basically, I am doing a logistic regression to estimate the length at which an individual matures at age t and t-1, then using an intermediate calculation to estimate the probability the individual is maturing for the first time (P_Mat <- (PrMat(t) - PrMat[t-1])/(1-PrMat(t-1)) in the code below). I'd then like to run a second logistic regression on P_Mat and length to determine the expected length at which an individual first matures. My attempted code is below:               model{          # Calculate Age-Specific Length at maturity with logistic regression (Ogive)     #Priors     for (i in 1:3) {    # Cohort 1-3         for (j in 1:4)  {   # Age 1-4             alpha[i,j] ~ dnorm(0, 0.0001)   # Intercept             beta[i,j] ~ dnorm(0, 0.0001)    # Slope         }     }          #Age and cohort-specific logistic regression centered for better convergence     for (i in 1:N) {             logit(p[i]) <- alpha[COHORT[i],AGE[i]] + beta[COHORT[i],AGE[i]] *                                      (LENGTH[i]-mean(LENGTH[]))             MAT[i] ~ dbern(p.bound[i])  # Bernoulli distribution for maturity              p.bound[i] <- max(0, min(1,p[i]))   # Draw from dist limited to 0-1     }          ##Intermediate calculations to derive P_Mat          for (i in 1:3) {          for (j in 2:4)  {              Prev_alpha[i,j] <- alpha[i, j-1]             Prev_beta[i,j] <- beta[i, j-1]         }     }     #Dummy variable for age-1 fish (not concerned about their PMRNs)     for (i in 1:3){         Prev_alpha[i,1] <- 1         Prev_beta[i,1] <- 1     }          #Estimate probability each fish matured at current age     for (i in 1:N){         PrMat[i] <- 1/(1+exp(-(alpha[COHORT[i],AGE[i]] + beta[COHORT[i],AGE[i]]*                     (LENGTH[i]-mean(LENGTH[])))))         PrevMat[i] <- 1/(1+exp(-(Prev_alpha[COHORT[i],AGE[i]] +                        Prev_beta[COHORT[i],AGE[i]](PrevLen[i]-mean(LENGTH[])))))         P_Mat[i] <- (PrMat[i] - PrevMat[i])/(1-PrevMat[i])     }          ######Begin model specification for second regression to estimate PMRN     #Priors     for (i in 1:3){         for (j in 1:4){             M_alpha[i,j]~dnorm(0, 0.0001)             M_beta[i,j]~dnorm(0,0.0001)         }     }          #Second regression of derived P_Mat data     for (i in 1:N) {             logit(p[i]) <- M_alpha[COHORT[i],AGE[i]] + M_beta[COHORT[i],AGE[i]                                (LENGTH[i]-mean(LENGTH[]))             P_Mat[i] ~ dbern(p.bound[i])  # Bernoulli distribution for maturity             p.bound[i] <- max(0, min(1,p[i]))   # Draw from dist limited to 0-1     }          }     Supplied data is AGE, LENGTH, PrevLen (estimated length of each fish        at previous age), MAT (1 or 0), and COHORT.      So, is this possible to do? If so, how? I was thinking wrapping the first regression and calculation of P_Mat in a datablock in JAGS could potentially solve the problem, but I've tried that to no avail. I've also read about cut and dsum functions, but I'm still learning the BUGS language and its possibilities, so any thoughts would be extremely helpful. Thanks!


title: Sum of lagged coefficients (VAR)

text: I have come upon a paper where a VAR analysis is performed. Using 3 endogenous variables and some exogenous (control) variables, the results of the VAR analysis are shown in a table. For the endogenous variables four lags are selected according to different selection criteria. However, in the estimation results table "the coefficients reported are the sum of the four lags; tests of their joint significance (F-stats) are reported in parentheses." My questions: 1) Can you sum (positive and negative) coefficients from a VAR analysis, and why? 2) How to sum and test for their joint significance in EVIEWS, if I wish to perform this operation myself? Thanks a lot!


title: Generalization of Newton-Raphson method for maximization of log-likelihood?

text: Suppose $l(b; y)$ is the log-likelihood function we want to maximize, i.e. we want to solve $U = l' = 0$ for the ML estimation of $b$. Using Newton-Raphson method, we have \begin{align*} b^{(m)} = b^{(m-1)} - \frac{U}{U'}, \end{align*} This can be again approximated by: \begin{align*} b^{(m)} = b^{(m-1)} - \frac{U}{E(U')}, \end{align*} If the random variable $y$ is from a distribution belonging to the exponential family, then we have $M = Var(U) = -E(U')$, thus: \begin{align*} b^{(m)} = b^{(m-1)} + \frac{U}{M}. \end{align*} If $\mathbf{b}$ is a vector of several parameters, and $\mathbf{U}$ is vector of derivatives where $U_j = \frac{dl}{db_j}$, then it's said the above formula can be generalized to \begin{align*} \mathbf{b}^{(m)} = \mathbf{b}^{(m-1)} + \mathbf{M}^{-1} \mathbf{U}, \end{align*} Where $\mathbf{M}$ is the variance-covariance matrix of $\mathbf{U}$ I am wondering how this generalization can be derived, book recommendations are also very appreciated.


title: How to visualize an enormous sparse contingency table?

text: I have two variables: Drug Name (DN) and corresponding Adverse Events (AE), which stand in a many-to-many relation. There are 33,556 drug names and 9,516 adverse events. The sample size is about 5.8 million observations. I want to study and understand the association/relationship between DN and AE. I’m thinking about a way of visualizing this set in R since it's better to look at pictures. I'm not sure how to go about it...


title: Generalized estimating equations output in SPSS

text: I am hoping to confirm that I have a suitable way to analyse the different proportions of people who are categorized as left lateralised on the one hand, or bilateral/right lateralised on the other in two different tasks. I cannot use an ordinary logistic regression (or chi square test) as the conditions are repeated measures. I have used the Generalized Estimating Equations option in SPSS to allow for the within subjects individual intercepts to vary as for repeated measures, but am wondering how best to interpret the output to show that the proportion of those in each category differs between the two tasks. Is it correct that I need to re-estimate the model but remove _condition_ from my fixed effects and then compare measures of model fit in each estimation for a significant difference (e.g.using the AIC)? Or should I just stick to looking for a significant Wald statistic and leave it at that? My final question is how to interpret the pairwise comparisons from this type of analysis. This shows a significant difference between my two conditions, but as I understand it to be a log odds value I'm not quite sure where to go with it. I do apologise if I haven't made my problems clear, I am quite new to this.


title: Kernel density estimation with an Epanechnikov kernel in MATLAB

text: I am currently trying to learn how to estimate the kernel density using the Epanechnikov kernel in MATLAB, and I am currently having problems with my code. What I'm doing is that the data that I am simulating comes from a mixture of normals. When I tried to estimate it using a Gaussian kernel, the code worked. However, when I used an Epanechnikov kernel, the graph I am obtaining is wrong. The code is attached here:               % Kernel density estimation using mixtures of normally distributed random     % variables.          % Set parameters for data generating process.     n=1000; % Number of observations.     data1 = rand(n,1); % Returns data with 1000 observations.          % Initializing the zero matrix of data from the mixture.     data2=zeros(n,1);          % % Generate the data. (mixture of normal random variables)     data2 = (data1 < ones(n,1)/3).*normrnd(2,1,n,1) + ...         (data1 > ones(n,1)/3).*normrnd(-2,1,n,1);          % Generating the density (using an Epanechnikov kernel.)     % Setting the parameters.     s=std(data2); % standard deviation of the data     h=2.34; % bandwidth parameter          % Evaluate the kernel at all x's in the domain.     x=linspace(-1,1,1000);     p=ones(n,1);     xi=x*p; % matrix for the x's where we evaluate the kernel.     data2i=data2*ones(1,n); % matrix for each of the data points in the kernel.     u=(xi-data2i)/h; % matrix of u's.     absu=abs(u); % absolute value of u.     I=(absu<=1);     f=(.75/h)*(1-u.^2).*I;          plot(x,f);      May I ask what I might be doing wrong?


title: Dealing with a categorical variable that can take multiple levels simultaneously

text: I recently posted a question with many parts and I'd like to focus in on just one issue that I didn't emphasize in the original post. My data is a list of records, each one representing an educational seminar event. I have a continuous variable that represents the revenue brought in by each seminar, which is the response variable in my regression. I also have a number of categorical variables which are acting as factors/IVs. One of those categorical factors is the speaker hosting the event. The trouble is that sometimes more than one speaker hosts a particular event. To date, all our speakers have been drawn from a pool of 154. Most of the time, just one speaker is used, but in about 10% of the data points, two, three or even four speakers were used. Currently, this is represented in my data with slashes ("Speaker One / Speaker Two / Speaker Three"). I've written a Python script that can find the average revenue on a given date interval for seminars which take a given level of a categorical variable (for example, it could return the average revenue for all seminars in 2008 that for which Speaker One was the host)...my script can read the multiple speaker format fine, reading names on opposite sides of a " / " as separate speakers. Unfortunately, R doesn't seem to be able to do anything like that...I've run a multiple regression on my data and obviously it treats "Speaker One", "Speaker Two" and "Speaker One / Speaker Two" as three different speakers. My multiple r-squared value is less than 0.5 so I'm hoping that resolving this issue would help improve the model...how best to proceed?


title: How to assess effect of intervention in one state versus another using annual case fatality rate?

text: I am a beginner in statistics with just basic knowledge. I have these data: cases, deaths and CFR (Case Fatality Rate-deaths per 100 cases) of a disease for 17 years (1994-2010) from 2 neighbouring states where people can walk across the states freely. This is a population based cohort study. Data are available from 1994. The treatment protocol was started in 1996. State 1 and state 2 implemented the same treatment but state 2 could not implement the treatment perfectly due to local administrative problems. The death rate fell in 1 state and continued to be high in the other state. Because it would be unethical to subject patients to do a case-control study, I want to analyse the following available data to see if there is a significant difference in the death rates and risk ratios of these 2 states because of poor implementation of the treatment guidelines in state 2 from 1996 to 2010. St. refers to State               Year   St.1 Cases  St.1 Deaths St.1 CFR St.2 Cases St.2 Deaths  State2 CFR  Risk Ratio     1994    1836        383         20.86     583        121         20.75        0.99     1995    1246        257         20.63    1126        227         20.16        0.98     1996    1450        263         18.14     896        179         19.98        1.10     1997    2953        407         13.78     351         76         21.65        1.57     1998    1161        149         12.83     1061       195         18.55        1.43     1999    2924        434         14.84     1371       275         20.06        1.35     2000    1729        169          9.77     1170       253         21.62        2.21     2001   1888      275         14.57     1005      199          19.80        1.36     2002    919     178          19.37      604      133         22.02         1.14     2003    865        142         16.42      1124      237         21.09         1.28     2004    1543       131          8.49      1030      228         22.14         2.61     2005    2887       336         11.64      6061     1500         24.75         2.13     2006    1484       108         7.28       2320     528          22.76         3.13     2007    1592       75          4.71       3024     645          21.33         4.53     2008    1920       53          2.76       3012     537          17.83         6.46     2009    1477       40          2.71       3073     556          18.09         6.68     2010    1534      26           1.69       3540     494          13.95         8.23      Kindly advise me what is the best way to go.   I have basic knowledge of using SPSS V19 and Comprehensive Meta-analysis V2.


title: When should I add a time-trend to a regression?

text: I've looked for this answer around the web with no luck so far. I'm mostly interested in how time trends apply to cross-section and panel data. Thanks


title: Trend or no trend?

text: I have a quarterly time series and test for stationarity with an augmented Dickey-Fuller test using R.               adf.test(myseries)     # returns     # Dickey-Fuller = -3.9828, Lag order = 4, p-value = 0.01272     # alternative hypothesis: stationary       so the H0 is rejected. I tried to validate this intuitively and regressed the same series on a linear trend.               x<- 104:1     fit.1<-lm(myseries~x)     summary(fit.1)     #returns     # x      0.024  1.31e-05 ***      Even though a simple linear model is not so appropriate here and the intercept is large (around 80), there seems to be a slight downwards trend over time, which is in line with my thoughts after looking at the initial data. So do I get the adf.test wrong or is the trend just to small to be discovered? Besides I used               plot(stl(myseries,"per"))      and ended up with a graph which sidebars suggested that trend and remainder were the main components driving the data, while seasonal influence was negligible. I saw that `stl()` uses Local Polynomial Regression Filtering and got a rough idea how that works (still I wonder why smoothed trends of Hadley's ggplot2 package looked that different even though it uses the same method by default). So summing up I got: - adf finding no evidence for a trend - a slight downwards trends "detected" by eyeballing and the naive approach - loess decomposition stating that the trend has strong influence (by the relation of its bars in the plot) So what can I learned from this? Probably I do have a terminology problem here, because the former two seem to address time trends while the latter address some other trend I cannot fully grasp yet. Maybe my question is just: Can you help me to understand the trend extracted by loess? And how is it related to smoothed / filtered stuff like HP-Filter or Kalman Smoothing (if there is a relationship and similarity does not only occur in my case)?


title: What is the proper ratio of mean squares for a two-level nested ANOVA?

text: I would like to analyse the following type of experiment using a two-level fully nested ANOVA. We have mice of two genotypes (Group factor) and for each of the genotypes we take a certain number of mice (sub-groups) and measure several identical samples per each mouse. Thus it looks like a classic two- level nested design with Genotype and Mice as fixed factors. I am interested in the main group effect of the genotype on the measurement outcome. I used _Statistica_ where F-value is calculated by the ratio of $\text{MS}_\text{group}/\text{MS}_\text{residual}$. However, I found several references on the web suggesting that it is better to calculate F-value from the ratio of $\text{MS}_\text{group}/\text{MS}_\text{sub-group}$. For example, on John McDonald's web page: > In a two-level nested anova, there are two F statistics, one for subgroups > (Fsubgroup) and one for groups (Fgroup). The subgroup F-statistic is found > by dividing the among-subgroup mean square, MSsubgroup (the average variance > of subgroup means within each group) by the within-subgroup mean square, > MSwithin (the average variation among individual measurements within each > subgroup). The group F-statistic is found by dividing the among-group mean > square, MSgroup (the variation among group means) by MSsubgroup. The P-value > is then calculated for the F-statistic at each level... On the same web page, he comments that in Rweb: > Fgroup is calculated as MSgroup/MSwithin, which is not a good idea if > Fsubgroup is significant. In my analysis, if there was a significant Group (Genotype) effect calculated as $\text{MS}_\text{group}/\text{MS}_\text{residual}$, I didn't observe a significant effect of sub-groups (Mice). Does it mean that I can use $\text{MS}_\text{group}/\text{MS}_\text{residual}$ ratio for estimating the F-value?


title: How to find the variance for a mean response in a multiple linear regression model

text: The question is the following: >   1. Our regression model can be written as $y_i = \beta^Tx_i + \epsilon_i, > 1 \leq i \leq n$. Find the $100(1- \alpha)\%$ confidence interval for the > mean response $\beta^Tx$ at a given vector $x$ of regressor observations. > >   2. Now suppose a future observation y is taken at a regressor $x$. Find a > $100(1-\alpha)\% $ prediction interval for $y$. > > Let:   * $k$ denote the # of regressors    * $n$ denotes the sample size    * $\sigma^2$ is ${\rm Var}(\epsilon|X)$ where $X$ is the $nxk$ matrix containing all the regressor values.   * $x_i$ is a $nx1$ vector of all the sampled values of the $i^{th}$ regressor  * * * So from my understanding the mean response is $E[y_i|x_i]$ where $1 \leq j \leq n$. And it is obvious that $E[y_i|x_i] = \beta^Tx$. However, I am unsure of how to get the variance needed to get the confidence interval for #1. Using $y_i \pm t_{\alpha/2, n-k-1}\sigma_{y_i}$. I used the following method as an attempt: $${\rm Var}(\beta^Tx_i|x_i) = x^T_i {\rm Var}(\beta|x_i) x_i$$ I took ${\rm Var}(\beta|x_i) = \sigma^2 (x^Tx)^{-1}$, which is the value ${\rm Var}(\beta_{OLS}|X)$, where $X$ is the $nxk$ matrix containing all the regressor values. So I got the answer: $$y_i \pm t_{\alpha/2, n-k}\sqrt{x_i^T\sigma^2(x^Tx)^{-1}x_i}$$ which is little more than a horrible guess since for a mean response only $x_i$ is given. I tried doing it this way too:   $${\rm Var}(\Sigma_{j=1}^k \beta_jx_{ij}|x_i)$$ but the expansion was too ugly to bear... * * * For #2, I got the following answer: $$ \hat y_{OLS} \pm t_{\alpha/2, n-k}\sqrt{\frac{e^Te}{n-k}} $$ which I am afraid another guess... where I assumed it was just an estimate using Ordinary Least Square method. Please point to a method of solving this variance problem. Please explain your rational!


title: How to draw a slope on a log-log scale in gnuplot

text: I have a log-log scale plot by gnuplot. With following ranges,               set xrange [1:1000]     set yrange [1:1000]      How do I draw slope 5/3 slope, ie, a slope that intersects y-axis at 5 and x-axis at 3? Many thanks!


title: Books for learning non parametric Bayesian model

text: Having studied parametric Bayesian statistics during the two last years, I plan to begin to self-study non parametric Bayesian model during this summer and look for recommendations. I would like the book(s) to cover both the theoritical aspects and the practical ones (with implementation, examples of models ...) Thanks for your recommendations.


title: Do multiple observers increase the probability of detecting an event?

text: I will describe my question with an example: A machine manufactures toys and places them on a conveyor belt. A toys expert overlooking the conveyor belt observes the toys. If a toy is faulty he can successfully detect it with probability $P_{d}$. It is possible to have more experts overlooking the process. If I have one expert and he detects a faulty toy, then I am $P_d$ sure that the toy is faulty (right?). Intuitively I understand that if more experts detect the toy as faulty then probability of the toy being faulty should increase (right?). What happens if 2 experts report the toy as fault and the other 2 report is as not faulty? What is the probability of the toy being faulty then? How can you work these questions out formally with probability theory?


title: Maximum likelihood estimation for custom equation in R

text: For this non-linear equation $$y = a\left(1-\exp\left(-\left(x/b\right)\right)\right)^c$$ by using least square estimated values for the parameters are obtained and               x <- c(3, 33, 146, 227, 342, 351, 353, 444, 556, 571, 709, 759, 836,      860, 968, 1056, 1726, 1846, 1872, 1986, 2311, 2366, 2608, 2676,      3098, 3278, 3288, 4434, 5034, 5049, 5085, 5089, 5089, 5097, 5324,      5389, 5565, 5623, 6080, 6380, 6477, 6740, 7192, 7447, 7644, 7837,      7843, 7922, 8738, 10089, 10237, 10258, 10491, 10625, 10982, 11175,      11411, 11442, 11811, 12559, 12559, 12791, 13121, 13486, 14708,      15251, 15261, 15277, 15806, 16185, 16229, 16358, 17168, 17458,      17758, 18287, 18568, 18728, 19556, 20567, 21012, 21308, 23063,      24127, 25910, 26770, 27753, 28460, 28493, 29361, 30085, 32408,      35338, 36799, 37642, 37654, 37915, 39715, 40580, 42015, 42045,      42188, 42296, 42296, 45406, 46653, 47596, 48296, 49171, 49416,      50145, 52042, 52489, 52875, 53321, 53443, 54433, 55381, 56463,      56485, 56560, 57042, 62551, 62651, 62661, 63732, 64103, 64893,      71043, 74364, 75409, 76057, 81542, 82702, 84566, 88682)          y <- c(1:136)     df <- data.frame(x,y)          fit <- nls(y ~ a*(1-exp(-x/b))^c, data=df, start = list( a=100,b=1000,c=0.5),            algorithm="port",lower=list(a=100,b=100,c=0.5),upper=list(a=200,b=10000,c=2))      **How to do this using Maximum Likelihood Estimation (MLE) for this non-linear equation?**


title: Probabilistic Max Pooling Convolutional RBM and ReLU unit?

text: In Honglak Lee's definition of a Probabilistic Max Pooling Convolutional RBM, the probabilities are formulated as: $$ P(h_{i,j}|v) = \frac{exp(I(h_{i,j}))}{1 + \sum\limits_{(i',j')\in B_\alpha}{exp(I(h_{i',j'}))}} \\ P(p_{\alpha}|v) = \frac{1}{1 + \sum\limits_{(i',j')\in B_\alpha}{exp(I(h_{i',j'}))}} $$ In this definition both the hidden and pooling units are binary. I was wondering if it could be adapted to use Rectified Linear Unit (ReLU) especially in the hidden layer since they have been proved more powerful ?


title: Applying the Akaike Information Criterion to Data

text: I have some variables that I would like to run regressions on, to create a model, but I am unsure about how to actually AIC (or the BIC). Unfortunately I have not yet taken a mathematical statistics course, so if someone could direct me to a nontechnical resource, or a resource that can explain how to use the methods without a formal background in statistics, this would be very much appreciated. Note: Though I have not taken a mathematical statistics course, my mathematics background is fairly strong (Calculus, Linear Algebra, Diff Eq, Probability, Elementary Analysis). So technical resources that don't presuppose a statistical background are not out of the question.


title: What data mining/machine learning approach to use for a scoring model?

text: Suppose I have a large data set with lots of features(attributes). And I'm tasked to build some kind of scoring model to rank certain objects with all these features. How do I go about doing this? From my understanding so far, I like to think of this as a supervised learning problem. But the problem is there is NO labeled classes (or at least it's not apparent). How can I rank order these objects? The closest thing I can think of is credit scores, but in credit scoring models, one supposedly has labeled classes as to who historically was good and bad. Should I invent/create some metric based on the list of attributes and use them as labeled cases? Like if attribute$_1 > x$ and attribute$_2< y$ etc., then it's considered "good" I believe they want a numerical ranking (i.e., scoring all the objects have numerical scores assigned to the objects like credit scores). If that's the case, then do I even need machine learning/data mining? Can't I just rank it by these attributes once they agree what the ordering means?


title: Covariances between forecasts from different time-series models

text: I'm trying to compute the variance of an average of forecasts. But I'm not sure on how to get the covariances required to compute the variances. Here is the situation: Three ARMA(p,q) models for a variable are estimated and using each model forecasts for 12 months ahead are produced. In all three models the residuals are normally distributed. An average forecast for each month is produced using the forecasts from ARMA(p,q) models. I would like to compute the variance for each month's average forecast in order to be able to provide the confidence intervals. The variance of sum of correlated variables are given by: $$Var(\sum_{i=1}^{n}{X_i})=\sum_{i=1}^{n}\sum_{j=1}^{n}Cov(X_i,X_j)$$ $$= \sum_{i=1}^{n}Var(X_i)+2\sum_{1 \le i}\sum_{1\le j}Cov(X_i,X_j)$$ Thus, the variance of an average forecast can be obtained by: $$Var(\sum_{i=1}^{3}\frac{X_i}{3})= \frac{\sum_{i=1}^{3}Var(X_i)+2\sum_{1 \le i}\sum_{1\le j}Cov(X_i,X_j)}{3^2}$$ Eviews spill out the variances for each forecast, but I'm not sure on how to get the covariances between the forecasts for each period. Any help would be appreciate! Thanks!


title: Cox model with Lasso for PH test and P value

text: I am using the `penalized` package in r for fitting a Cox model with a lasso penalty. Does the Cox model + lasso need to test proportionality? How do you get the p-value for the nonzero coefficients in the model? Thanks for any response.


title: Discrimination and calibration of Cox model

text: I have been working on fitting Cox model for prediction by using the rms package. I want to measure model calibration and discrimination. Discrimination was measured by using `rms::validate()`; `Dxy` can be transferred to @Frank Harrell's $c$ index. But in this way, I cannot get 95% CI for the $c$ index. How can I do this in R? And by the way, what value should be in $c$ index to present the model's well? Calibration was done using `rms::calibrate()`, but I cannot get the calibration plot which presented concordance of predicted and observed events in Cox model. How can I do this calibration plot of Cox model in R or SAS.


title: How to conceptualize error in a regression model?

text: I am attending a data analysis class and some of my well-rooted ideas are being shaken. Namely, the idea that the error (epsilon), as well as any other sort of variance, applies only (so I thought) to a group (a sample or whole population). Now, we're being taught that one of the regression assumptions is that the variance is "the same for all individuals". This is somehow shocking to me. I always thought that it was the variance in Y accross all values of X that was assumed to be constant. I had a chat with the prof, who told me that when we do a regression, we assume our model to be true. And I think that's the tricky part. To me, the error term (epsilon) always meant something like "whatever elements we don't know and that might affect our outcome variable, plus some measurement error". In the way the class is taught, there's no such thing as "other stuff"; our model is assumed to be true and complete. This means that all residual variation has to be thought of as a product of measurement error (thus, measuring an individual 20 times would be expected to produce the same variance as measuring 20 individuals one time). I feel something's wrong somewhere, I'd like to have some expert opinion on this... Is there some room for interpretation as to what the error term is, conceptually speaking?


title: Testing for difference between incidence rates in R, spatial rdd

text: I am an economics student just starting out with R, and while I'm beginning to be somewhat comfortable with it, I also realize I need to strongly brush up my basic statistics. I have individual level admin data (in the form of a rotating panel) documenting a reform which was implemented as a pilot project in some regions but not in others, and I want to analyze the effects of the reform in a spatial regression discontinuity design. I have already extended the data with distances and travel time which I acquired via the google maps api. To get a first impression, I have compiled incidence rates in treated and non- treated regions for different time frames, and different distances to the nearest resp. non-treated/treated location. This is an excerpt of the results table:                    sample base  ref treated distance  all    recip         irr     1    pp.1yr 2003 2004       0    1e+04  602        5 0.008305648     2    pp.1yr 2003 2004       0    3e+04 6357       39 0.006134969     3    pp.1yr 2003 2004       0    5e+04 8528       57 0.006683865     4    pp.1yr 2003 2004       0    1e+05 9272       62 0.006686799     5    pp.1yr 2003 2004       1    1e+04  435        4 0.009195402     6    pp.1yr 2003 2004       1    3e+04 2438       16 0.006562756     7    pp.1yr 2003 2004       1    5e+04 3456       22 0.006365741     8    pp.1yr 2003 2004       1    1e+05 6360       45 0.007075472     9    pp.2yr 2002 2004       0    1e+04  245        2 0.008163265     10   pp.2yr 2002 2004       0    3e+04 2693       25 0.009283327     11   pp.2yr 2002 2004       0    5e+04 3699       36 0.009732360     12   pp.2yr 2002 2004       0    1e+05 4084       37 0.009059745     13   pp.2yr 2002 2004       1    1e+04  187        1 0.005347594     14   pp.2yr 2002 2004       1    3e+04  983       11 0.011190234     15   pp.2yr 2002 2004       1    5e+04 1400       16 0.011428571     16   pp.2yr 2002 2004       1    1e+05 2660       35 0.013157895      (I can also reshape everything to wide.) My primary question is simply practical: **How do I test in R whether the corresponding rates for treated/non-treated differ significantly, i.e. whether rates in treated regions are significantly smaller than in untreated regions?** I was looking at `t.test()`, but I guess `prop.test()` is appropriate here. From eyeballing it doesn't look like there is any effect in the direction I had expected. Any further info on how to proceed and hints at problems I might encounter or apparent misunderstandings would be nice. One thing I thought of is that the number of observations is quite small. The main reason besides the rather small target group is the anonymisation of locations in the data, which leaves me with a lot of useless observations because I cannot locate them. There is a distant possibility of this restricion being relaxed so that I can use more data to increase the sample size (esp. in border regions), but I am rather sceptical. I also remembered to look at pre-/post- differences before treatment and after implementation nationwide. Any help is greatly appreciated.


title: Justification of one-tailed hypothesis testing

text: I understand two-tailed hypothesis testing. You have $H_0 : \theta = \theta_0$ (vs. $H_1 = \neg H_0 : \theta \ne \theta_0$). The $p$-value is the probability that $\theta$ generates data at least as extreme as what was observed. I don't understand one-tailed hypothesis testing. Here, $H_0 : \theta\le\theta_0$ (vs. $H_1 = \neg H_0 : \theta > \theta_0$). The definition of p-value shouldn't have changed from above: it should still be the probability that $\theta$ generates data at least as extreme as what was observed. But we don't _know_ $\theta$, only that it's upper-bounded by $\theta_0$. So instead, I see texts telling us to assume that $\theta = \theta_0$ (not $\theta \le \theta_0$ as per $H_0$) and calculate the probability that this generates data at least as extreme as what was observed, but only on one end. This seems to have nothing to do with the hypotheses, technically. Now, I understand that this is frequentist hypothesis testing, and that frequentists place no priors on their $\theta$s. But shouldn't that just mean the hypotheses are then impossible to accept or reject, rather than shoehorning the above calculation into the picture?


title: Sphering data with SVD components of covariance matrix

text: The Elements of Statistical Learning says on page 113: > Sphere the data with respect to the common covariance estimates > $\hat{\Sigma}$: $X^* \leftarrow D^{-1/2}U^TX$ where $\hat{\Sigma} = UDU^T$. > The common covariance estimate of $X^*$ will now be the identity. Can someone help me understand why $X^* \leftarrow D^{-1/2}U^TX$ spheres the data?


title: How to generate random variates from random variables with known density?

text: Given a random number generator to produce random variables with a probability density function $f(x)$, how to generate random variables with probability density function $g(x)?$


title: The unit information prior and its BIC approximation

text: Just moments ago, I asked this question because I've been reading Wagenmakers 2007. I now have a better understanding of what a unit information prior is and can push my knowledge further with two (by which I mean 3) questions:   1. The unit information prior has been considered by some to be too conservative. Are there other criticisms unique to the unit information prior that are not covered by criticisms of "objective priors" as a family? And, does this prior compare favourably to similar "objective priors"?   2. It was raised in the previous question that there is some _discomfort_ using the BIC approximation. I would appreciate an elaboration as to why this may be.


title: model with truncated count response using negative binomial regression model

text: 1. what is the clear model for truncated negative binomial regression model with brief explanation about the parameters?   2. what is the clear difference between truncated negative binomial regression model and zero-truncated negative binomial regression model?


title: wrong reported Total Sum of Squares in time fixed effects with plm (twoways)

text: I was sent here from Stackoverflow because this is more a statistics question. I hope I am in the right place now! The summary command of a plm regression with the (factor="twoways") argument reports the same coefficients as a plm regression with manual time dummies (and the default factor="individual") but a much lower total sum of squares (and therefore a much lower R2). Why is this the case? **Example:**               library(plm)          individual <- c(rep("AT",5),rep("BE",5),rep("DE",5))     time <- c(rep(2006:2010,3))     set.seed(123)     y <- c(rep(rnorm(15, sd=2)))     x <- c(rep(rnorm(15, sd=3)))          Test <- data.frame(individual=individual, time=time, x=x, y=y)     year.f = factor(Test$time)     dummies = model.matrix(~year.f)     Test <- cbind(Test,dummies)     remove(dummies,year.f)          fe_manual <- plm(y~ x+year.f2007+year.f2008+year.f2009+     year.f2010,data=Test,model="within")     summary(fe_manual)          fe_twoways <- plm(y~ x, data=Test,model="within",effect="twoways")     summary(fe_twoways)      The problem doesn't look very bad here (we go from total sum of squares of 38.7 in the manual model to 30.1 in the twoways model), but in my real sample, this problem gets me from an R2 of 42% to one of 3%. The problem gets much bigger with more data. edit: You can see that the twoways model calculates time fixed effects by running:               fixef(fe_twoways)      Strangely, these coefficients are again different from the manual model.


title: Model selection with glmer

text: I need to do a backward elimination on my glmer model, but I have no idea to do that ? It seem that a lot of thing do no work with glmer ... thank you


title: Ranking based on multiple parameters

text: I am xyz organization. I have recruit entry level engineers from 100 colleges. I have their performance, attrition & promotion data for last 3 yrs. I have 4 parameters performance, promotion, attrition and volume intake. I want to rank these colleges based on the above parameters using any statistical model or tool.


title: Finding repeating patterns in sensor data using Java

text: I'm working with a dataset that has accelerometer (tri axis) readings, magnetometer (also tri axis) readings, readings for colour etc. All the values are numerical. This data was collected by wearable camera technology. I am trying to find repeating patterns in this data, but I'm new to data mining and not sure what the best approach would be. I have read about clustering algorithms, and nearest neighbour also. However as I say, I'm very new to this and not sure what approach is best. If anyone has any suggestions on where to start it would be a big help. Thanks!


title: Precision of an incomplete classifier

text: Given a testing set of nodes which can be either +,-, or 0, I use an incomplete classifier which allows me to predict if a node is +, -, 0, not +, not -, or not 0, and sometimes it cannot predict anything. How are precision, recall, and accuracy defined for such a classifier?


title: Outlier detection in binary classification

text: I have a question about outlier detection in my system. I’m designing a system (in Matlab) that optimize both features and parameters of a classification method (like mlp) together with optimization algorithms. My input data are financial statement and some these statements have strange and unusual values. (For example someone use only round values like 10000, 100000 in every field or used 1 or 2 in some fields. So I want use outlier detection in my system. Besides it is an option in my system to use dimensionally reduction (like FDA – my system is binary classification) for preprocessing. My critical questions is here:   1. Now I’m using k-means clustering for finding outliers. I cluster every binary class [0 or 1] and remove samples with highest distance from center (for example in every cluster 5% of data). Is this a good approach for outlier detection in my binary classification system? Which outlier detection techniques do you suggest in this case?   2. Where should I put outlier detection function? Before data normalization (mapmaxmin or mapstd), after it? Before or after dimensionally reduction? Or finally because my system find best inputs in among 21 inputs of my system (financial ratios) in every iteration for chosen inputs I must outlier detection? (I think maybe these outliers leaving negative impacts on normalization or outlier detection)   3. After training my system and find best system structure I want use it in real with new dataset and find class of inputs [0 and 1]. So in this phase, I must use outlier detection for this data? If yes, where should I put outlier detection in this phase? Before normalization, after it? Before dimensionally reduction or after it? Thanks


title: Finding best parameters of SVM in matlab

text: I’m designing a system (using Matlab) that I can optimize parameters of a support vector machine (SVM) with genetic algorithm, harmony search and another optimization algorithms to find the best structure of SVM for a specific data. My problem is binary classification with 0 and 1 output and I normalize data (mapmaxmin o mapstd) before insert it to system. Besides it in some cases I use dimension reduction (for example FDA) to reduce my features. For this normalized data I must set the boundary of searching space in optimization algorithm. This is my SVM function:               svmstruct=svmtrain(TrainInputs,TrainTargets,...     'boxconstraint',Penalty,...     'kernel_function','rbf','method','QP',...     'rbf_sigma',Sigma,...     'autoscale','false');      I optimize only 'boxconstraint' and ‘rbf sigma’. For boxconstraint, my algorithm is searching in [0.001 400] and for sigma the searching space is same [0.001 400]. IS this searching boundaries is suitable for my problem or I must change these boundaries? Otherwise, I set ‘autoscale’ to ‘false’. Which one is better in my problem? ‘false’ or ‘true’ ? I set kernel function to rbf. is this a good approach for this problem? Thanks


title: Multicollinearity if independent variables sum to one?

text: I am intending to explain the variation in dependent variable Y through a number of ratios that, combined, sum to 1 in each period `t` Hence, $$ Y_t = \alpha + \beta_i^I * X_t + \epsilon_t $$ Where $X_t$ is a vector of ratios that sum to 1 (say 1/3, 1/6 and 3/6 in t = 1) and $\beta_i^I$ a vector of its beta coefficients. In actually, the ratios reflect the composition of a fund, that changes over time. Is there multicollinearity if I include all the ratios? That is, should one of the ratios be excluded from the model (such as with the dummy variable trap) to avoid multicollinearity? Since the ratios change disproportionally to each other but proportionally to the sum of 1, I'm not sure if this should be the case.


title: Functions for regression diagnostics on mer objects in R

text: I'm particularly interested in plotting residuals against fitted values, and residuals against predictors. Often times I need to make boxplots of the residuals conditional on predictors. I'd be interested in functions that include other/multiple diagnostics as well, e.g. ceres plots, qqplots. The car package has some nice diagnostic functions but they are only for objects of class lm and I often use lmer.


title: Robust standard errors for panel data vs robust estimation for panel data

text: It seems to me that rather than using non robust estimation methods with robust standard errors it would be better to use robust estimation from the outset. I wonder what other people think.


title: Is there a notion of the product of confidence bands?

text: Given two confidence bands, it seems natural to consider the confidence band of the product of these two distributions. Does this have a name? Are there any reference that might be helpful regarding this notion? Edit: Assume that the errors of the first, say, regression curve are not correlated with the errors of the second.


title: ARIMA model fitting and forecasting implementation

text: I am new to time series modelling and I am trying to build a simple time series model using ARIMA methodology and forecast. I could write an R program to do the same, but I am more interested in writing my own `arima(1,1,0)` and `predict()` functions that R provides in Java / C++ / Python. While I can go about reading the theory behind ARIMA modeling and then implementing it, I still face a learning curve there, due to time constraints, I would like to take an approach where I can see a time series data sample, and a write up which walks me through the steps to compute `arima()` and `predict()` so that I can implement an algorithm. If someone can point me to a sourcecode / psedocode / step-by-step example via data, it would really be helpful.


title: Modelling a skewed, 10-point Satisfaction variable

text: I am trying to replicate and hopefully improve on an analysis done in a study to find determinants of patient satisfaction after shoulder surgery. Satisfaction is heavily skewed (with over 60% of patients responding 10) The original study ran a multiple linear regression on satisfaction with a handful of predictors. The residuals are just as skewed as the satisfaction variable itself, and no transformations I've tried have helped, probably because the variable is discrete (could argue ordinal or numerical). My other thoughts have been to use ordinal regression on the full 10-level scale, or lump responses together into two or more groups for binary (say 10 vs not 10) or multinomial logistic regression. The former seems to suffer from decreased interpretability and the later seems to lose information that might be important. Am I missing any ideas/techniques? How do I choose between these various evils? Thanks in advance for any responses!


title: What test to asses the complication a treatment has with a small sample size(80) and only 5-10% complication rate?

text: I am assessing the benefit of using a therapeutic intervention in a population of 80 people split evenly into two groups. The benefit is to decrease a certain complication which only has a incidence rate of 5-10%. How can I prove non-inferiority? Which statistical test is best an and how can I calculate the sample size needed for power? I assume this is a pretty easy question, hope someone can help! Thanks


title: Large overlap between confidence intervals, although z test for difference was significant

text: I have calculated a z test to examine the difference between two proportions. I found differences to be significant. The formula I used was: $$z = (p_1-p_2)/SE$$ where $$SE = \sqrt{ p ( 1 - p ) ( \frac{1}{n_1} + \frac{1}{n_2} ) }$$ and $n_1$ is the sample size of sample 1, $n_2$ sample size of sample 2. I have also plotted a bar graph of my proportions with 95% confidence intervals around each of the two proportions based on this equation: For proportion 1:   $p_1 \pm 1.96SE$ $SE = \sqrt{p(p-1)/ n_1}$; $n_1$: size of the sample 1 For proportion 2: $p_2 \pm 1.96SE$ $SE = \sqrt{p(p-1)/ n_2}$; $n_2$: size of the sample 2 I would not expect the CIs to overlap much based on significant of z test. However, they overlap by more than 25%. Am I using the correct formula for the confidence intervals? If so why does there seem to be discrepancy between the results of the z test and the visual depiction of the difference


title: In R, what does a probability density function compute?

text: I understand that the probability density function, pdf, of a continuous random variable is the probability of the variable taking on a given value. I am also thought that for a continuous random variable, the probability of it taking a specific value is always 0 i.e. P(X=x) = 0. Therefore, in R, I am wondering what does > dnorm(0) compute? I know it means the probability of obtaining a 0 from a normal distribution (which is one of the distributions a continuous random variable can take) with mean 0 and standard deviation 1. I had expected the output to be 0, since I am asking for the probability of a continuous random variable taking on a specific value. However, the output is 0.3989. What does 0.3989 represent and why isn't the output 0?


title: Violation Proportionality Cox model - Repeat analysis?

text: I'm new to survival analysis, but I've been reading some papers and books and I got a nice model. However, one of the variables (Sit) does not met the proportionality assumption for the Cox model. Nevertheless, this got me thinking and it makes sense, since that variable should have a time-varying effect. For instance, in this case, I expect that shortly after treatment, the risk for relapse is higher for some individuals while as time passes, the time for all individuals is similar. Note (to myself) that time-varying effect is different than time- varying covariate. It got me quite some time to fully understand this. Now, my **question** is: I know how to include an interaction between Sit and time (I'm using R). But when I add the interaction some of the variables that were previously significant are not anymore. And it got me thinking, **should I repeat the process to select variables?** Perhaps, some of the variables previously removed are now significant. Or should I interpret the model as it is, and ignore the variables that are no longer significant? Another **question** is: The proportionality of the Cox model is violated. I have decided to introduce an interaction. **So, now I should not do the residuals analysis?** Since the proportionality is no longer in order. **How should I test for the fitness of the model, then?** I'm really confused about this and it would be great to have some guidance and tips about the best procedure, from now on. Thank you very much!


title: Difference in predicted value using two different methods

text: Take these two vectors:               points <- c(44, 36, 33, 33, 29, 28, 28, 22, 21, 20, 18, 15, 15, 15, 12, 12, 12, 11, 10, 10,  8,  8,  8,  8,  7,  7,  7,  6,  6,  6,  6,  6,  5,  5,  4,  3, 3, 2,  2,  2,  2,  1,  1,  1,  1,  1,  1,  1)          hours <- c(137.000000,  58.450000,  92.250000,  94.750000,  80.000000, 100.000000,  33.750000,  48.500000,  70.500000,  76.500000, 43.250000,  90.750000,  33.000000,  13.750000,  16.250000,  27.250000,  30.830000,  12.000000,  68.000000,  34.000000, 13.250000,  32.250000,   9.750000,   4.083333,  18.000000,  23.750000,   4.666667,  31.750000,  11.750000,   4.850000, 10.866667,   4.166667,  14.000000,   6.166667,   2.000000,   7.750000,  11.100000,   4.750000,   1.750000,   0.250000, 1.000000,  12.000000,  13.000000,   1.000000,   0.250000,   5.250000,   5.333333,   2.166667)      I made this model:               lm(points ~ hours)          # Call:     # lm(formula = points ~ hours)          # Coefficients:     # (Intercept)        hours       #     2.9530       0.2823      The model suggests that if I do 1 hour, I get 3.2353 points. I have also calculated mean points per hour:               mean(points/hours)          # 0.839875      The ratio suggests that if I do 1 hour, I get 0.839875 points. Why is there such a big difference in the number of points per hour predicted by each method? Or am I interpreting something wrong?


title: Upper Bound of Measured Data

text: HW question I am having trouble on: ![enter image description here](http://i.stack.imgur.com/qY8nh.jpg) My attempt at solving it: ![enter image description here](http://i.stack.imgur.com/uh9zh.jpg) I thought part (a) was pretty straight forward until I worked it out and got the same value of K for K_upperbound which I assume is incorrect. Also how would I apply the least squares method to this problem?


title: How to use SPSS for this statistic question? Please help me :(

text: I am in the basic level of statistic... I have some trouble in how to do SPSS for this statistic question. So my topic is if music with or without lyric would make a impact on reading comprehension task. I have two same music but with and without lyric. And also I have two reading (A & B) comprehension that is about the same length and it contain 5 questions each. But there will be order effect on which music is presented first and which reading comprehension question presented first. So my teacher suggested to do 4 different group with   1. with lyric (ask to solve A reading comprehension first) then without lyric (reading B)   2. With lyric (ask to solve B reading comprehension first) then without lyric (reading A)   3. Without lyric (ask to solve A reading comprehension first) then with lyric (reading B)   4. Without lyric (ask to solve B reading comprehension first) then with lyric (reading A) I kind of know what i am doing but not clear how to make a excel file (I mean how to organize this stat)and put it in to SPSS (still not familiar with this programm) Since my volunteer will participate in both condition it would be a paired independent sample t test? Our independent variable is music played at two different conditions, with and without lyric. The dependent variable is the performance score, which will be measured by the score of reading comprehension task.


title: Given the density of X and a sample size of 4 with yielded observations 1, 2.5, 3.5, 4 calculate the maximum likelihood estimator

text: I'm currently just learning about Maximum likelihood estimators and I'm a bit stuck on how I should answer this question. ![a busy cat](http://i.imgur.com/lzlS1Xu.png) So far I know that for the function to hold, the observation 4 should be <= 2*theta. and that theta>=2 as a result right? But then I'm not sure how to use it to find the estimator. I set the Likelihood function as the following ![a busy cat](http://i.imgur.com/sNpXisP.png) I know I'm supposed to differentiate the log likelihood functions with respect to theta but I get -1/theta=0. Is there no mle estimator or am I missing something?


title: Does R predict or forecast on an arimax model?

text: I have some data which i am trying to work on. I am pretty new to R though, but i love R. First, I fit an `arima` model to this data and used the `detectIO` function in R to detect a single influential outlier (IO). I then incorporated this IO into my model and then developed an `arimax` model, now with the IO. I first used this model:               Model2 = arimax(mydata, order=c(0,2,1),                     seasonal=list(order=c(0,0,1), period=12), io=c(48))      Then later, used the model below after seeing a similar code in a book:               Model4 = arima(mydata, order = c(0, 2, 1),                    seasonal=list(order = c(0, 0, 1), period = 12),                     xtransf=data.frame(I48=1*(seq(mydata)==48), I48=1*(seq(mydata)==48)),                     transfer=list(c(0,0), c(1,0)))      I then run the `predict` function on these two models and I keep getting this error:               > predict(Model4, n.ahead=24)          Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),  :       attempt to set an attribute on NULL      I don't know what I am doing wrong. I don't know if I am using the right functions and codes? Does there really exist a `predict` function on `arimax`? I would also appreciate if I get any help on how to set-up the `xtransf` with some real data values.


title: What is the meaning of a large p-value?

text: I understand that the $p$-value is the conditional probability of observing the test statistic or something more extreme given that the null hypothesis is true. I have read the great explanation by @user28 in this post: What is the meaning of p values and t values in statistical tests? However, do large $p$-values say anything? Does a larger $p$-value lend greater support to the null hypothesis? If I set rejection region to be $<0.05$, then does it make a difference if I get $p$-value $0.06$ or $0.99$? (After all, $0.05$ is arbitrary, and $0.06$ is so close to being rejected that if I arbitrarily set $0.05$ as $0.1$ instead, the null hypothesis would have been rejected.) Can one make any statistical use of a non-rejecting $p$-value?


title: What is the justification for unsupervised discretization of continuous variables?

text: A number of sources suggest that there are many negative consequences of the discretization (categorization) of continuous variables prior to statistical analysis (sample of references [1]-[4] below). Conversely [5] suggests that some machine learning techniques are known to produce better results when continuous variables are discretized (also noting that supervised discretization methods perform better). I am curious if there are any widely accepted benefits or justifications for this practice from a statistical perspective? In particular, would there be any justification for discretizing continuous variables within a GLM analysis?    * * * _[1] Royston P, Altman DG, Sauerbrei W. Dichotomizing continuous predictors in multiple regression: a bad idea. Stat Med 2006;25:127-41_ _[2] Brunner J, Austin PC. Inflation of Type I error rate in multiple regression when independent variables are measured with error. The Canadian Journal of Statistics 2009; 37(1):33-46_ _[3] Irwin JR, McClelland GH. Negative consequences of dichotomizing continuous predictor variables. Journal of Marketing Research 2003; 40:366–371._ _[4] Harrell Jr FE. Problems caused by categorizing continuous variables. http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/CatContinuous, 2004. Accessed on 6.9.2004_ _[5] Kotsiantis, S.; Kanellopoulos, D. "Discretization Techniques: A recent survey". GESTS International Transactions on Computer Science and Engineering 32(1):47–58._


title: When sample size is small and distribution is not normal should I use Wilcoxon or T-test or something else?

text: I have two sets of data and need to compare mean times for a variety of situations (i.e. time of day, day of the week, different areas, etc...). Some of the situations have small samples sizes (n=8) and need to be compared against a larger sample size (n=260). Neither data set is normally distributed (both are usually positively skewed). The larger data set has a standard deviation that is approximately 3 times the size of the mean. The p-value for the T-test is 0.000003. The p-value for the Wilcoxon is 0.61. It seems to me that given the nature of the data Wilcoxon may be my best bet, but the p-value for the T-test seems more believable. Am I approaching this correctly?


title: Forecasting using multiple regression

text: I have data in the form given below, and I want to perform forecasting using multiple regression. I found definition of multiple regression from this link: http://otexts.com/fpp/5/1/ . I have these questions: **(1)** Can we forecast value of Y (the dependent variable) based on the given values of the independent variables (X1 and X2), as in below table? or it is mandatory to have previous values of Y along with X1 and X2 to forecast the future value of Y, **(2)** I want to forecast the value of Y at T8, by providing two inputs, i.e., value of Y_at_T8 = forecast(X1=22, X2=33), how to proceed on that? , **(3)** Is it mandatory for multiple regression that Time must be in fixed intervals, can it increase in varying intervals? **Time | X1 | X2   ** T1 | 23 | 33   T2 | 43 | 22   T3 | 42 | 44   T4 | 33 | 32   T5 | 43 | 11   T6 | 55 | 21   T7 | 23 | 42   Regards!


title: What statistics should I use for evaluating the accuracy of predictions?

text: I have two variables representing 1) players' predicted fantasy football points and 2) players' actual fantasy football points scored. What statistics are best for assessing the accuracy of the predictions in terms of absolute agreement (not just relative rank-order accuracy) in relation to the actual values? I'd like to compare the accuracy of various sources of projections. The predicted and actual values are both continuous. I'm considering statistics like: r-squared, Harrell's c-index, and ICC (absolute agreement). Which statistic(s) would you recommend and why? Thanks in advance!


title: Rescale factor scores from factor analysis to latent metric in R

text: I'm calculating a factor analysis of several variables in R. I want to determine each case's value on the latent variable. When I run the factor analysis, I receive factor scores. The factor scores do not have the same metric as the raw data or the latent variable. How can I rescale the factor scores to correspond to the metric of the latent variable (same mean, sd, and range) to determine each case's value on the latent variable? Here's a small example:               v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)     v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)     v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)     v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)     v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)     v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)     m1 <- cbind(v1,v2,v3,v4,v5,v6)          factanal(~v1+v2+v3+v4+v5+v6, factors = 1, scores = "Bartlett")$scores      Thanks in advance!


title: using random forest for missing data imputation in categorical variables ( in R)

text: I have following type of **associated** data. The following example step to generate associated variable. **p** number of variables and **n** is number of observations.               p = 500     n = 200     mat <- matrix(NA, ncol = 500, nrow = 200)          for (i in 1:p){     if(i ==1){          fs <- sample (c("AA", "AB", "AB", "BB"), n, replace=TRUE)           mat[,i] <- fs           fs1 <- fs           }      rechr <- sample(1:n, 1)     fs1[rechr] <- sample (c("AA", "AB", "AB", "BB"), 1)     mat[,i] <- fs1      }          dim(mat)     mat[1:10,1:20]      The above data is complete, but in real dataset I have missing values (which are randomly distributed within each variables). I would like to predict them using random forest (or any other appropriate algorithms). Let's randomly put 10% missing values (the approximate number of variables) to the above data.               rowind <- sample(1:n, 20)     colind <- sample(1:p, 20)          for (i in 1:length(rowind)){                     mat[rowind[i],colind[i]] <- NA     }     mat[1:20,1:10]      How can do this and **what issues** I need to think up in doing so - considering **categorical variables** , number of observations and variables ? **Edit:** As all variables are correlated and none are response (but can serve as response) while prediction. for example mat[,1] missing value prediction can be based on rest of variables, mat[,2] prediction can be based on prediction of rest of variables. So the objective here is prediction all missing values in all variables (all of them are categorical). **Edit 2** I am interested in multiple imputation of this data so that I can perform other analysis that requires complete dataset without missing value.


title: Robust MCMC estimator of marginal likelihood?

text: I'm trying to compute the marginal likelihood for a statistical model by Monte Carlo methods: $$f(x) = \int f(x\mid\theta) \pi(\theta)\, d\theta$$ The likelihood is well behaved - smooth, log-concave - but high-dimensional. I've tried importance sampling, but the results are wonky and depend highly on the proposal I'm using. I briefly considered doing Hamiltonian Monte Carlo to compute posterior samples assuming a uniform prior over $\theta$ and taking the harmonic mean, until I saw this. Lesson learned, the harmonic mean can have infinite variance. Is there an alternative MCMC estimator that is nearly as simple, but has a well-behaved variance?


title: Why don't we square the reliability correlation when calculating reliable variance?

text: This question relates to situations where we apply some test to the same people at Time 1 and Time 2, and then calculate the correlation between scores at Time 1 and Time 2. I am reading the book Principles & Applications of Assessment in Counseling, which states here on p.52 that: > Therefore, if the reliability coefficient using test-retest were .80, using > classical test theory we would interpret it by saying that 80% of the > variance is true to observed variance and 20% is error to observed variance. > In reliability, we don't square the correlation coefficient or use the > coefficient of determination. Instead, we simply use the reliability > coefficient itself to evaluate the degree of measurement error. I've seen it suggested in many places that it's a common mistake to square the measurement error in this scenario. However, I do not understand why it is a mistake or precisely what is meant by 'reliable variance'. Whiston, S. C. (2009). Principles & applications of assessment in counseling. Belmont, CA : Brooks/Cole, Cengage Learning.


title: Assessing Cannibalization, intervention of new mobile app on monthly sales

text: I am a beginner in statistics and looking for suggestions from you all on the approach for one of my study. For my study, there is a company which sells products via its online website (lets call it Desktop Website). They have recently launched a mobile app in addition to the desktop website. Both online website & mobile app sell the same stuff. **What I am try to understand:** 1) Whether sales through mobile app is cannibalizing sales through online website 2) Whether introduction of mobile app has led to increase in the total sales ( through both online website & mobile app) **I have the following in my hand:** 1) Monthly sales data from 2011 Jan (to date ) for our Desktop Website 2) Monthly sales data from 2012 Sep (to date ) for our mobile app. The mobile app was launched only on Sept 2012 3) Page views on mobile & Desktop website for the aforementined respective dates. Pakages: SPSS, Excel. **General observations from the data:** 1) The sales (for both mob & online) shows **annual seasonal effects** ( sales more during holiday months) 2) The sales (for both mob & online) shows a **positive linear trend**. **What I did:** 1) To understand intervention effect of the mobile apps on the desktop sales, I regressed the Desktop Website sales data with 11 seasonal dummies and an intervention dummy ( 0 for months before Sept 2012 and 1 otherwise) _Results:_ The R square was 67% (Adjusted R square =53%). But Durbin-Watson showed .789 which meant autocorrelation exists on the residual errors. So even though the intervention_dummy coefficient had highest value and with p value almost 0%, I had to reject the model because of the auto correlation. 2) In the next model, I included the mobile sales figures too. i.e I regressed the Online sales data with Mobile sales data and 11 seasonal dummies ( Did not include intervention dummy because SPSS automatically removed it while running the analysis ) _Results_ : R square = 95%, Adjusted Rsquare= 91%. Durbin-Watson = 2.17 & the Regression Standardized Residual histogram shows almost a perfect normal distribution. The coefficient for Mobile Sales figure was low and corresponding pValue was at 40%!. **_Does this mean that mobile sales is not impacting the desktop sales ?_**. Can I safely conclude that cannibalization is not taking place? What do yo suggest for assessing the impact of mobile app on the overall sales? Preferably I am looking for something that can be done with help of statistical packages and less of heavy mathematical models if that's possible. But any suggestion is welcome as I think I have hit analysis paralysis! Thanks a bunch!


title: Predicting the year computers will be able to do 1 exaFLOPS using historical data (time series forcasting using prediction interval)

text: Put simply, I need to make a prediction interval (or confidence interval?) for the x value when y reaches a certain number given historical data for x versus y. Long Version: I have data on the speed of computers in FLOPS from the 1950s up until 2012 and need to make a prediction interval for the predicted year we should reach an exaflop (10^18 FLOPS). My problem is that prediction intervals require means and standard errors, and all I have is a single year that my initial regression predicted along with the data used to predict that value. Now, It's clear that that single predicted value is not sufficient to make any significant statistical inferences on. So, what I think I need to do is use my historical data to construct a confidence interval for the year a computer should reach an exaflop so that I can run a hypothesis test, where the Null Hypothesis is that the year we should reach an exaflop is 2020 and the Alternative Hypothesis is that the year we should reach an exaflop is >2020. The sort of end result I'm thinking of is something along the lines of being able to say "I am 95% Confident that the true year a computer should be able to do an exaflop is between 2020 and 2035", in addition to the statements for the hypothesis test. My issue now is that what I THINK I need to do is construct many regression lines on the data to get a mean predicted year for the exaflop and the accompanying standard error so that I can do the prediction interval. Is this possible? If so how do I do it? If there exists a better, easier way to predict the year of an exaflop given my historical data please let me know what it is and how to do it. Help is greatly appreciated. EDIT: After researching time series forecasting, I've become aware that what I'm trying to do is construct a prediction interval on a time series. I haven't found a reasonable and accepted way to carry this out yet. Help is still needed and appreciated.


title: Purpose of Scheffé's confidence intervals

text: Scheffé's method, when first written about by Henry Scheffé, was described as serving the purpose of identifying statistically significant contrasts in multiple-comparison's ANOVA problems. But here in this forum I am told in answers and comments under this question that that should not even be done: "A test of contrasts will only have the correct distribution under the null when the contrasts are specified before you see the data." As stated, that is of course nonsense. It depends on what test is used. If you use at test designed only for testing whether a particular contrast differs significantly from zero, and you find that there are _some_ contrasts that are significant, then you're making a mistake. But that's what I take Scheffé's method to be intended to remedy. And in a comment: "And in general it's unwise to even consider testing patterns in data based solely on the pattern in the data." What is Scheffé's method _for_ if not exactly that? It seems to me that some people have been instructed that testing a contrast suggested by the data is not valid if one uses the same test one would use if that particular contrast had been identified in advance (which is obviously right) and leapt from there to the conclusion that there are no _correct_ methods of testing contrasts suggested by the data. "The naive method is incorrect; therefore no method is correct." ---- my paraphrase Could those who want to defend the view I'm paraphrasing explain what they think the purpose of Scheffé's method is?


title: Deriving a "confidence relationship"?

text: I understand the basics of confidence intervals, the central limit theorem, etc, to be able to know things like given N samples of random variable, we're 68/95/99.7 percent sure the variable is within these two values. I'm wondering about how this generalizes to cases where we're trying to derive a relationship between two variables. Ie, if we're trying to measure how some variable Y depends on a variable X, one thing I guess we could do is set X to 1, then get a 95% confidence interval for X = 1. Then set X += 0.1 and do that again, and again again etc, to derive a 95% "confidence relationship", or something. This seems kind of sloppy though. Is there some other statistical mechanism designed for doing things like this (ie, finding a relationship between X and Y that we are 95% sure is the right relationship). edit: Some more info on what I'm looking for: Say I have a bunch of (x,y) pairs. I'm trying to express the y's as a function of x. x values might repeat. I also don't necessarily care about expressing the relationship algebraically. So, the naive thing to do would be to average over the ys for every possible x value we have data for, and plot that. (so if the dataset was {(1,2), (1,3), (2, 2)}, I'd plot the points {(1,2.5), (2,2)}). This totally disregards data about how much y data we got for each x value, though. It also seems to be ignoring the fact that we suspect the data to be correlated in some way. Like if our data points were {(1 + 1e-10, 4), (1, 17)}, my naive approach would plot {(1,17), (1 + 1e-10, 4)}. But considering how close 1 and 1e-10 are, we'd think that the values 4 and 17 should be close to each other. So maybe something like {(1, 10.5), (1 + 1e-10, 10.5)} would be more accurate. Or, in another sense, if our dataset was {(1,16), (1+1e-10,16.2}), the fact that we got two values so close to 16 both times should seem to "inspire our confidence". But since the x values are technically different we won't notice this.


title: Identifying a minimum time on web page as proxy for successful search outcome

text: Given a correlation between two sequences of numbers, I need to find to what point in one sequence the correlation is strongest. Here's why: I'm tracking large numbers of pages in a knowledge base, and flagging certain pages for followup. For example, if a page is among the most frequently viewed but also has a high ratio of negative ratings, it makes sense to prioritize that for further investigation and improvement. One of the more specific things I'd like to target is pages that have a high number of unrelated searches. It's no good if many people land on a page, but it's not actually what they were looking for. In that case, perhaps the page title and first paragraph or meta description need to be improved to cut out misleading vocabulary and more clearly indicate what the page is about. However, it's very time-consuming to look at search queries for hundreds of pages. I'd like to get a proxy metric to flag potential problem pages first, before diving into the search queries for those pages. I think I've found such a metric: average time on page. After manually grading a large number of pages based on the ratio of related keywords to unrelated keywords, I found that there was a positive correlation between low time on page and a high ratio of unrelated keywords. To be precise, an Excel `CORREL` formula came back with 0.46, which I'm led to believe is a reasonable indicator of a correlation. (By now, I'm sure it's obvious that I'm not a professional statistician. ;-) ) Is this a good start? If so, what kind of Excel formula could I run on the existing data to identify a minimum time on page, below which a page would be flagged for further followup? (I.E., a formula that would find the best cutoff point, above which the time on page didn't correlate so strongly with unrelated searches.)


title: Guessing data from statistics

text: If I have a number of statistics/information of a data set, e.g. mean, SD, skew, kurtosis, % data values within a certain range, etc, how might I go about trying to deduce a portion of the actual data. In my case, the data are not numbers, but rather subjects in a study, and therefore have more information (and thus associated constraints). I just don't know where to begin with this problem. Is there some branch of statistics that might be relevant to look into? Some have pointed me to sampling, but this seems quite general. Thanks!


title: Proof of density inequality

text: I was wondering if there is an easy way to find sufficient conditions for the following inequality to hold $$ \int f(x,y)^2 \:\mathrm{d}x \:\mathrm{d}y - \int f(x)^2 f(y)^2 \:\mathrm{d}x\:\mathrm{d}y \geq 0, $$ where $f$ is the density function. Do you have any suggestions?


title: Coefficient standard error of zero in quantile regression

text: I'm currently experimenting with quantile regression of a strongly right skewed outcome variable y on a 3-category exposure x (values 1,2,3). I wanted to model the .2, .5, and .8 quantile, using the interior point algorithm for estimation and bootstrapping (Markov chain marginal bootstrap method, n=1000) for confidence intervals. The software (SAS v9.2) gave no warning about convergence problems or the like. For the latter two, the results appear reasonable, however, for the .2 quantile, I sometimes get standard errors of 0 for the intercept and both of the two exposure coefficients (dummy coded, 2 vs. 1 and 3 vs. 1). Group sizes for subgroups defined by x are: 227 (x=1), 423 (x=2), and 90 (x=3). Like I said, y is strongly right skewed, with a range of 0-14 (mostly integer values, but decimal values also exist in the data). Proportion of subjects with y=0 is 31% for x=1, 18% for x=2, and 10% for x=3. I speculated that this might be the reason for the 0 standard errors, but I neither have enough experience with this particular method, nor enough mathematics background to really understand the algorithms. Therefore my questions: Is it normal to get stanard error of 0 with this method? And what does it imply for the interpretation of these results? Should I drop the analysis of the .2 quantile altogether? Does the validity of the other two quantile regressions (for .5 and .8) remain unaffected by this? Below is the SAS code I used, and the results of an exemplary analysis. I repeated the analysis 10 times and in 5 cases, standard error were 0 like in the example below, and in the remaining 5, they were between .01 and .10. Results for the .5 and .8 quantile were pretty stable. PROC QUANTREG DATA=... ALGORITHM=INTERIOR(TOLERANCE=1e-4) CI=RESAMPLING(nrep=1000) ORDER=DATA; CLASS x; MODEL y = x / QUANTILE = .2 .5 .8; RUN; I'd be highly grateful for any comments!               Quantile Parameter      x value   Estimate    StdErr       LowerCL      UpperCL          0.2      Intercept                0.0000      0.0000       0.0000       0.0000     0.2      x              3         1.0000      0.0000       1.0000       1.0000     0.2      x              2         1.0000      0.0000       1.0000       1.0000     0.5      Intercept                1.0000      0.0805       0.8419       1.1581     0.5      x              3         3.0000      0.4168       2.1817       3.8183     0.5      x              2         1.0000      0.0828       0.8375       1.1625     0.8      Intercept                3.0000      0.2288       2.5508       3.4492     0.8      x              3         4.2019      0.5462       3.1296       5.2742     0.8      x              2         1.0000      0.2444       0.5203       1.4797


title: What test should be used to tell if two linear regression lines are significantly different?

text: I have two lines in a x-y plot:   1. observed migration distance (y) and predicted migration distance (x) $$Y=0.95X+0.31,$$ where 0.95 and 0.31 are the slope coefficient and the intercept (p-value = 0.001)   2. 1:1 line (slope: 1 and intercept: 0) The null hypothesis is $$y=ax+b$$ $$a_1=1=a_2\quad and\quad b_1=b_2$$ Please kindly advise what test to use to test the above hypothesis.


title: How to analyse qualitative interviews coded on a 1 to 4 ordinal scale?

text: I'm a student helping a professor with some research interviews done in QSR Nvivo. Basically the interviews are coded in 4 outcome categories:   1. solely positive   2. generally positive with a few considerations   3. generally positive with some considerations   4. generally negative  On top of this coding are a lot of subcodings, that I would like to crosscheck with the The sample consists of ~60 interviews. Only about six interviews were in category 4 (negative) . **Questions:**   * How do I run some stats on this ordinal scale (1-4) according to the subcategories?   * Is it a valid scale to use? (I hear that comparing 2 v1 is pretty unusual, but it has been done, but I haven't seen 3v1.)   * Can I get anything statistically significant based on a between group 20v20v14v6 sample?


title: varimax rotation, Direct Oblimin and factor analysis

text: What is the difference between varimax rotation and direct oblimin in factor analysis? Also, I am confused the relationship among principal component analysis, varimax rotation and exploratory factor analysis in theory and SPSS? Are they have some connections? Thank you.


title: How to statistically compare groups for multiple density plots?

text: Is there a statistical method to compare these density plots other than ANOVA (MANOVA)? I would like to compare the densities among year within each plot and report which of those distributions are "significantly different". The reason is that I am more interested in the spread and where on the distribution of the dataset rather than solely on the mean and variance. Thanks Here is the code and plot               ggplot(NMPSCFAM, aes(Length, colour= Year)) +        geom_density(aes(y = ..count..), alpha=0.3) + xlim(0, 700) +         geom_density(aes(colour = Year)) +       xlab("") +       ylab(expression(paste("Mean density ( ", m^2, ")", sep = ""))) +     facet_wrap( ~ Family + Sector2, ncol=3, scales = "free")      ![enter image description here](http://i.stack.imgur.com/vI6pI.png) * * * Edit: Ok so no answers so far. If I did these density plots in a histogram would there be a way to do a Chi-square test on each of the plots? Could any one help me out on how to do this in R?


title: what test to use tell significant difference of two quadratic regression curves

text: I have two curves of quadratic regression models in a x-y plot   1. $y=a_1\text{x}^2+b_1\text{x}+c_1$ $Y= -0.51\text{x}^2-0.88\text{x}+3.21$, $R^2: 0.12$, coefficient of the quadratic term: p-value = 0.001   2. $y=a_2\text{x}^2+b_2\text{x}+c_2$ $Y= -0.17\text{x}^2-0.13\text{x}+3.41$, $R^2: 0.99$, coefficient of the quadratic term: p-value = 0.001 I want to know if curve 1 and 2 are significantly different. The null hypothesis is $a1=a2$ and $b1=b2$ and $c1=c2$ Please kindly advise what test to use to test the hypothesis


title: Degrees of freedom of J-Test distribution

text: This book states, on page page 256 ( the GMM section) that the J-test(for over-identifying restrictions) is of the following form $G_n'J_n^{-1}G_n \approx \chi^2(m-p) $, where $G_n=G_n(\theta)=1/n \cdot \sum^n_{i=1}g_i(\theta)$ which are the sample versions of DGP original $m$ distinct moment conditions $E(g_i(\theta_0))=0$. The author says that this result comes from using $\frac{1}{\sqrt{n}}G_n\rightarrow_dN(0,J_0)$, and $J_0\approx \frac{1}{n} J_n$. I'm trying to deduce the degrees of freedom of the chi-squared distribution. For that I'm using the following result: Let $y_1,...,y_n\sim N(0,1)$, then if $y=(y_1,...,y_n)'$ and $A$ a symmetric and idempotent $n\times n$ matrix, we have that $y'A y\sim \chi^2(rank(A)) $. So, we know that $\frac{1}{\sqrt{n}}J_0^{-1/2}G_n(\theta_0)\rightarrow_dN(0,I)$, thus, $G_n'J_n^{-1}G_n \approx \left(\frac{1}{\sqrt{n}}J_0^{-1/2}G_n(\theta_0)\right)'\cdot I\cdot \left(\frac{1}{\sqrt{n}}J_0^{-1/2}G_n(\theta_0)\right)\rightarrow_d \chi^2(rank(I))$, where $I$ is an identity matrix of dim $m\times m$. So, I'm doing something wrong. What is it? Any help would be appreciated.


title: Random slope in a model without random intercept

text: I am looking for a model to fit longitudinal data and opted for GLMM with a logit link. I tried to add random effects but when both a random slope and random intercept are in the model, this does not converge. Then I tried to include only one of them and both models converged but I am wondering if I can have a model with random slope without a random intercept. Does this conflicts with the hierarchy of the model? Or, since all subject start at the same point in my data (it's about roses and at day 0, they are all fresh), is it allowed? The AIC value does drop when the random effect is a random intercept only but I am not sure if I can compare the AIC values. Thank you for any help!


title: Interpreting Omitted Variables

text: I do not fully understand how to interpret the difference between two statistical models where they only differ based on whether a certain variable is included on the right hand side. If the results do not change much, then do we say that that this omitted variable has little effect on the result? If the coefficient on the variable of interest changes significantly in magnitude and is no longer significant, we say that the variable of interest on the right hand side is not robust to that omitted variable. That omitted variable effects the relationship between the variable of interest on the right hand side and the dependent variable. What is the estimate becomes more statistically significant when this variable is included? It appears that including these variables makes the results more precisely estimated. How can this be? What does that mean about the relationship between the dependent variable and the main variable of interest on the right hand side if the relationship is more precisely estimated only with the additional control variable?


title: Is there a R command for testing the difference in coefficients of two linear regressions?

text: I am looking for a R command to test the difference of two linear regressoon betas. Lets say I have data $x_1, x_2...x_{n＋1}$. $\beta_1$ is obtained from regressing $x_1$ to $x_n$ onto $1$ to $n$. $\beta_2$ is obtained from regressing $x_2$ to $x_{n＋1}$ onto $1$ to $n$. Is there a way in `R` to test whether $\beta_1$ and $\beta_2$ are statistically different? Edit: Here is a clearer problem statement: I changed the notation for data from $x$ to $z$... That's it. Should be very clear now... Thanks! Data observations: $z_1, z_2, ..., z_{n+1}$               y1 = z_1,z_2,.........  z_n      y2 = z_2, z_3,......... z_{n+1}          x1 = 1, ..., n     x2 = 1, ..., n          y1 = A1+ x1 * B1 + epsilon_1     y2 = A2 + x2 * B2 + epsilon_2      H0: B1 and B2 are statistically significally different...


title: Comparison of stock market volatilities

text: I was trying to compare the volatilities of the stock market index of China and America. But I realized that a direct comparison of historical data would not be fair. The Chinese market has a rule that restricts changes in stock prices to be within $\pm$ 10% everyday. No such rule exists in the US. I couldn't just eliminate days where the American market changed by more than 10%, because the next day the change would have been bigger if previous' day's change was restricted to be within 10%. So is there any fair way to modify the data for comparison?


title: How to combine predictions in ensemble

text: I am trying to learn more about how to build ensembles of predictions in R and coming to a roadblock, and am hoping one can offer guidance. I often read about people automatically identifying how they should weight each model through the use of OLS. How do people do this? Do you just insert your prediction from each model as a regressor in the model? E.g., Final_Prediction = b0 + b1*prediction_from_GBM + b2*prediction_from_SVM + bkxk + e and just fit the line above to combine your predictions? What about when you are predicting class membership.... do you get the probabilities from each model and fit them in a logistic model similarly? Any resources that are R specific or any thoughts / clarifications are greatly appreciated. I have not been able to find anything.


title: Cohen's d and multiple comparisons for 2/3-way ANOVA

text: I am conducting three-way ANOVA (A*B*C) with 2 levels each.   1) I found A*B interaction.   2) I moved to 2-way ANOVA (A*B) and found interaction again. I reported Eta- squared and equivalent Cohen's d.   3) I moved to 1-way ANOVA (a1,a2,b1,b2) with mushed cells and bonferonni posthoc showed that b2 was lower than the a1,a2, and b1. Question 1: In step 2, what does Cohen's d tell me?   Question 2: In step 3, did I need to divide my alpha=0.05 by what?   Question 3: Can I report Cohen's d between b2-a1 and b2-a2 and b2-b1? thanks


title: How do I fit a simple filter with R?

text: I have two timeseries, `y[t]` and `x[t]`. A domain expert tells me that `y[t]` should be a linear combination of past values of `x[t]` up to a certain horizon, i.e.               y[t] = a[0]*x[t] + a[1]*x[t-1] + ... + a[n]*x[t-n]      How do I estimate the `a[i]` coefficients in R? Indeed, what is the proper name for this kind of model?


title: Can structural equation models be used to derive clinical formulae?

text: Structural equation models (sem) are used to model latent variables. Renal function is a latent variable measured by serum creatinine levels (with measurement errors) expressed by many different clinical formulae derived from linear regression models. Inulin clearance is a gold standard for renal function. However, hundreds of articles had studied the relative "accuracies" of these formulae as compared to surrogate "gold standards". Can sem be used to derive a clinical formula from serum creatinine and inulin clearance for the estimation of renal function?


title: When correlation coefficient's value rises, error rises as well. Is this normal behaviour?

text: So here's a question. lets say we have three datasets A, B, C, D(continuous, that take values within the data range [0, 5]) where C and D are actually the predicted values of A and B accordingly, taken from a model use. Both pairs (A and C, B and D) are significantly correlated and Spearman's Rho for pair (A,C) is 0.2 where Rho for pair (B,D) is 0.4. When computing the MAE though, the pair with biggest Rho also appears to have the maximum error. Shouldn't it be the other way round? I mean, the biggest the correlation coefficient, the smallest the error since we're talking about positively correlated data? Thank you in advance :) A.


title: Is there a generalization of trimean to $n$-th order (central) moments?

text: I think trimean is the cat's meow. Is there a generalization of this idea to $n$-th order (central) moments? Basically I live in a world where the pain of outliers vastly exceeds the pain of inefficiency. The reference calls this "macrostatistics" although I've never heard this term until researching this question.


title: Odds of X occurrences in a row given Y trials (A coin flip problem)

text: I 'm looking to find a general formula for the following problem: The events:   * Let each event have two possible outcomes (e.g. success and failure)   * The events are independent Scenario:     * Let the number of trials be Y   * Let the number of successes be N Here's the tough part: Over Y trials, what is the probability of there being X successes in a row given that the total number of successes is N?


title: What is the difference between serial correlation and having a unit root?

text: I may be mixing up my time series and non time series concepts, but what is the difference between a regression model that exhibits serial correlation and a model that exhibits a unit root? In addition, why is it that you can use a Durbin-Watson test to test for serial correlation, but must use a Dickey-Fuller test for unit roots? (My textbook says this is because the Durbun Watson Test cannot be used in models that include lags in the independent variables.)


title: Maximizing Cronbach alpha to select the number of variables that define a given construct

text: I can't seem to find this anywhere on the internet or in my books. I'm putting together a construct for an underlying variable. Is it acceptable to have a construct with two variables if that produces a bigger Cronbach alpha than the construct with 3 variables in a professional statistical analysis? Using 2 variables seems a bit on the low side, but it produces better results.


title: How to remove stopwords with Russian documents?

text: I have a bit of a problem with the tm R package for cleaning text documents. Here is how my code looks like:               library("tm")          # import text files in corpus text – ok     c_txt <- Corpus((DirSource(directory = ".", pattern ="txt",                      encoding = "UTF-8")), readerControl = list(language = "rus"))          # convert to Lower Case – ok     txt_cl <- tm_map(c_txt, tolower)          # remove Stopwords - not working     txt_cl <- tm_map(txt_t, removeWords, stopwords("russian"))      I am inspecting the file `russian.dat` which is located in the subdirectory `tm\stopwords`. There I can notice strange characters -- not Russian stopwords, although it is UTF-8 encoding. Where is my mistake? Thanks in advance.