title: Expected mean square of random effect in CRD model

text: Random effect in CRD model: \begin{align} y_{ij} &= \mu + \tau_{i} + \varepsilon_{ij} \\\ \tau_{i} &\sim \mathcal N(0,\sigma_{\tau}^2) \\\ \varepsilon_{ij} &\sim \mathcal N(0,\sigma^2) \end{align} I want to prove $$ E(MStr)=\sigma^2+n\sigma_{\tau}^2 $$ The following is what I have done so far: $\text{Let }i = 1, \ldots, p,\quad\quad j = 1, \ldots,n.$ \begin{align} E(\overline{y_{i.}}) &= E\bigg(\frac{\sum_{j=1}^{n}y_{ij}}{n}\bigg) = \frac{\sum_{j=1}^{n}E(\mu+\tau_{i}+\varepsilon_{ij})}{n} = \mu+\frac{\sum_{j=1}^{n}\tau_{i}}{n} \\\ \ \\\ &= \mu+\overline{\tau_{i}} \\\ \ \\\ E(\overline{y_{..}}) &= E\bigg(\frac{\sum_{i=1}^{p}\sum_{j=1}^{n}y_{ij}}{np}\bigg) = \frac{\sum_{i=1}^{p}\sum_{j=1}^{n}E(\mu+\tau_{i}+\varepsilon_{ij})}{np} \\\ \ \\\ &= \mu+\frac{\sum_{i=1}^{p}\overline{\tau_{i}}}{p} \\\ \ \\\ E(MStr) &= \sigma^{2}+\frac{n\sum_{i=1}^{p}(E(\overline{y_{i.}})-E(\overline{y_{..}}))^2}{p-1} = \sigma^{2}+\frac{n\sum_{i=1}^{p}\big(\overline{\tau_{i}}-\frac{\sum_{i=1}^{p}\overline{\tau_{i}}}{p}\big)^2}{p-1} \\\ \ \\\ &= \sigma^2+n\sigma_{\overline{\tau}}^2 \end{align} Could anyone help me find the mistakes I made?


title: What statistical tests for headache journal?

text: I track my pain levels in an online spreadsheet along with daily habits and trigger events. I want to test whether changes in my pain over time follow a trend (not concerned whether it is linear or not, just whether it is going up or down). I also want to test which daily events and triggers can be significantly associated with a change in pain level. It's been a while since I took statistics class and I'm confused by online flowcharts. Which statistical tests should I use and should my data be considered parametric or non-parametric? Thanks.


title: When "additive modeling" is better?

text: Given $y=f(x)+\epsilon$ where $x=(x_1,\dots,x_p)$, $f$ is highly non-linear and two different estimators:   1. $\hat{y}=\hat{M}(x)$   2. $\hat{y}=\hat M_1(x)+\hat M_2(x)$ where $M_1$ is a simple (biased) model and $M_2(x)=E[y-\hat M_1(x)|x]$. Can you tell me when the second estimator (two-step modelling) will give better predictions ? If $f$ is highly non-linear, are $M_1$ and $M_2$ easier to estimate than $M$ ? What about the bias/variance of the two estimators ? They have the same bias but what about the variance ? Is there a link with additive modeling ?


title: Need help trying to do PCA and then Gower Similarity

text: I have a dataset and i have reduced it using PCA, but once i reduce the dimension of the dataset, i lose all my uniqe columns in the data. Which i think is what should happen. But i want to run gower similarity on this and dont know how to run this using SPSS. I have looked at the macro online but cannot run it. Anyone has any idea on how to run similarity measures for data with different kind of data. I am new to this and LOST :(


title: Error message with MANOVA

text: I am having trouble with a MANOVA. I have an experiment with 6 response variable (say, A, B, C, D, E and F), and 3 explanatory variables (say R1, R2 and R3). I used `cbind()` (which I will call "y") to link my 6 response variables together, and then used this in my model as follows:               model <- manova(y ~ R1*R2*R3)      This is working without problem, I can obtain the results from my MANOVA, everything is fine. However, as my response variables are habitats of different sizes, it has been pointed out to me that I should correct my data to account for these differences in size. Once I do this, however, my MANOVA sends an error message which I do not understand :               Error in summary.manova(model) : residuals have rank 5 < 6      It seems that this message usually appears when a dataset has more columns than rows, which is not my case (I have about 9 columns and about 132 rows). Also, it is confusing me as my dataset after correction is very similar to the one before the correction (same amount of columns, same amount of rows, and in both cases all the values for my response variables are between 0 and 1). Any idea why I get this error message with the second model and not the first one?


title: Does it make sense to include combined averaged values as part of a crossvalidation round?

text: I have a large set of sequences S. A subset of these sequences are Recombination Signal Sequences (RSS). We have a subset of these sequences which we know are RSS. We are trying to predict which of the other sequences in S are also RSS. We try to show this prediction works by using a crossvalidation round. We divide the RSS into n disjoint subsets. We using (n-1)/n of the RSS as the training set, and (1/n) as the testing set. For each of the n values, we calculate a score (positive real number) for all the sequences in S. We rank the sequence in ascending order of score. The smaller the score, the more likely it is that the sequence is a RSS. For n=5, I graph the top 40 values for each of the 5 crossvalidation procedures in the round. This is shown in the figure below. My question is: the figure includes a sixth panel, which is an attempt to summarize the crossvalidation round in one plot. This averages the values over the different crossvalidation procedures. However, I was wondering if this really makes sense. Is there any sensible way to average/summarize these values with a suitable interpretation, or is this plot best left out? ![enter image description here](http://i.stack.imgur.com/gmuSE.png)


Q1: How can i convert Relative risk into Absolute risk? Q2: How to estimate absolute risk for an individual with the available sources of RR and OR? pl provide with example or related material thanks Raja


title: Help about a perceptron question

text: while studying for my Machine Learning exam, I encountered a problem that I cannot understand. In the problem, we have this perceptron, which 3 binary inputs (0 or 1) a,b,c with respective weights of 1,2,1 and a bias of 1.5. The question is what is the boolean logic function implemented by the binary perceptron and the answer is b + ac, (b OR (a AND c)), but I cannot get why this is the correct answer. My logic is that following the rule where the net sum has to be > threshold (1.5) in this case, then if the input a and c are on, then the perceptron should fire 1, which according to the answer is wrong since b has to be on. Can anyone help me understand this? Thanks.


title: Normalization Factor Wrong? (Bug?)

text: I'm new to PyMC and Bayesian stuff in general, so I started off with what I thought was a very simple toy problem. I generated some normally-distributed noise with a given mean and standard deviation, and then played with PyMC and looked at the posterior distribution of possible mu/sigma. It's possible that I just don't understand what I'm looking at, but shouldn't the standard deviation of the trace for mu converge on the standard error of the mean? It's currently off my a factor of sqrt(2). Simply: Parent normal distribution with parameters mu0, sigma0. Child population sampled from parent distribution, resulting in mu1, sigma1. SE = sigma1/sqrt(N) SD = sigma0/sqrt(N) Shouldn't SE~SD~mcmc.trace['mu'].std() for sufficiently large N and # of iterations? What I actually see is: SE~SD~mcmc.trace['mu'].std()/sqrt(2) This is using PyMC 2.2 and here is the code:               %pylab inline     import pymc as mc          true_noise_mu = 50     true_noise_sigma = 5     true_noise_datapoints = 2000          noise_field = np.random.normal(true_noise_mu, true_noise_sigma, true_noise_datapoints)               data_mean = noise_field.mean()     data_std = noise_field.std()               mu = mc.Normal( 'mu', data_mean, 1/(data_std/sqrt(n_datapoints))**2)          # just guessed at tau here, but it shouldn't impact the result (nor does it), #only the time to converge, yes?     sigma = mc.Normal( 'sigma', 1/data_std**2, 0.01)               observation = mc.Normal( "obs", mu, sigma, value=noise_field, observed=True)          mcmc = mc.MCMC(model)     mcmc.sample( 400000, 20000, 1 )          mu_samples = mcmc.trace('mu')[:]     sigma_samples = mcmc.trace('sigma')[:]          print "SEM estimated from sample population: %.3f" % (noise_field.std()/sqrt(true_noise_datapoints))     print "SEM calculated from parent population: %.3f " % (true_noise_sigma / sqrt(true_noise_datapoints))     print "SEM calculated from posterior distribution: %.3f" % (mu_samples.std())     print "Ratio of posterior: sample estimate: %.3f" % ((mu_samples.std())/(noise_field.std()/sqrt(true_noise_datapoints)))     print "1/sqrt(2): %.3f" % (1/sqrt(2))      Last Print:               SEM estimated from sample population: 0.677     SEM calculated from parent population: 0.671      SEM calculated from posterior distribution: 0.478     Ratio of posterior: sample estimate: 0.707     1/sqrt(2): 0.707      I'm probably just an idiot, but I would love to understand what I'm missing. =)


title: Type III sums of squares

text: I have a linear regression model with one categorical variable $A$ (male & female) and one continuous variable $B$. I set up contrasts codes in R with `options(contrasts=c("contr.sum","contr.poly"))`. And now I have Type III sums of squares for $A$, $B$, and their interaction (A:B) using `drop1(model, .~., test="F")`. **What I am stuck with is how sums of squares is calculated for $B$. I** think it is `sum((predicted y of the full model - predicted y of the reduced model)^2)`. The reduced model would look like `y~A+A:B`. But when I use `predict(y~A+A:B)`, R is returning predicted values that are the same as the full model predicted values. Therefore, the sums of squares would be 0. (For the sums of squares of $A$, I used a reduced model of `y~B+A:B`, which is the same as `y~A:B`.) Here is example code for randomly generated data:               A<-as.factor(rep(c("male","female"), each=5))     set.seed(1)     B<-runif(10)     set.seed(5)     y<-runif(10)          model<-lm(y~A+B+A:B)          options(contrasts = c("contr.sum","contr.poly"))          #type3 sums of squares     drop1(model, .~., test="F")     #or same result:     library(car)     Anova(lm(y~A+B+A:B),type="III")          #full model     predFull<-predict(model)          #Calculate sum of squares     #SS(A|B,AB)     predA<-predict(lm(y~B+A:B))     sum((predFull-predA)^2)           #SS(B|A,AB) (???)     predB<-predict(lm(y~A+A:B))     sum((predFull-predB)^2)      #Sums of squares should be 0.15075 (according to anova table)     #but calculated to be 2.5e-31          #SS(AB|A,B)     predAB<-predict(lm(y~A+B))     sum((predFull-predAB)^2)               #Anova Table (Type III tests)     #Response: y     #             Sum Sq Df F value Pr(>F)     #(Intercept) 0.16074  1  1.3598 0.2878     #A           0.00148  1  0.0125 0.9145     #B           0.15075  1  1.2753 0.3019     #A:B         0.01628  1  0.1377 0.7233     #Residuals   0.70926  6


title: Feature selection without target variable

text: Let's assume I have a NxD matrix X with the N rows being observations and the D columns being features. I would now like to know which are the most "interesting" features of this dataset. I.e. which features depend on each other, which are redundant etc. At the end, I would like to have a dataset of dimensionality k < D, because I could dismiss (D-k) features. My first idea was using PCA to get an approximation to the "intrinsic" dimensionality of my dataset. However, PCA will not directly tell me which features are the most interesting ones, it will only give me a number of principal components and their "strengths" (eigenvalues of the covariance matrix of X). So I thought about using a classical feature selection method like stepwise regression. However, stepwise regression requires a target vector y (since it is regression, of course) which I don't have. I only have the dataset X. **I only have basic machine learning skills, so I would like to know what is the appropriate method to select the most interesting features of my dataset X without having a target vector y.**


title: Multiple regression with "overlapping" categorical dummy variables

text: Using R, I want to run a linear regression to estimate the abnormal return on days with positive, negative and neutral news (CLASS). I'm a beginner in R, as well as in using regression models! First of all the data structure is as follows. CONTROLVAR just represents all the columns I use as control variables.               DATE <- c("1","2","3","4","5","6","7","1","2","3","4","5","6","7")     COMP <- c("A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B")     RET <- c(-2.0,1.1,3,1.4,-0.2, 0.6, 0.1, -0.21, -1.2, 0.9, 0.3, -0.1,0.3,-0.12)     CLASS <- c("positive", "negative", "aneutral", "positive", "positive", "negative", "aneutral", "positive", "negative", "negative", "positive", "aneutral", "aneutral", "aneutral")     SUBJECT.1 <- c("LITIGATION","LAYOFF","POLLUTION","CHEMICAL DISASTER","PRESS RELEASE","PEOPLE","EMISSIONS","ENERGY","WASTE MANAGEMENT","EMPLOYEES","SUBJECT11","SUBJECT12","SUBJECT13","SUBJECT14")     SUBJECT.2 <- c("POLLUTION","EMPLOYEES","NUCLEAR","FUELS","STOCK OPTION PLAN","EXECUTIVES","CO2","SOLAR","POLLUTION","EXECUTIVES","SUBJECT21","SUBJECT22","SUBJECT23","SUBJECT24")     SUBJECT.3 <- c("ENVIRONMENT","JOB REDUCTIONS","POWER PLANTS","POLLUTION","EMPLOYEES","FRAUD","CLIMATE CHANGE","SUSTAINABILITY","HAZARDOUS WASTE","BONUS PAY","SUBJECT31","SUBJECT32","SUBJECT33","SUBJECT34")     CONTROLVAR <- c("11","13","13","14","13","14","12","11","13","13","14","13","14","12")          df <- data.frame(DATE, COMP, RET, CLASS, SUBJECT.1, SUBJECT.2, SUBJECT.3, CONTROLVAR, stringsAsFactors=F)          df          #    DATE COMP   RET    CLASS         SUBJECT.1         SUBJECT.2       SUBJECT.3 CONTROLVAR     # 1     1    A -2.00 positive        LITIGATION         POLLUTION     ENVIRONMENT         11     # 2     2    A  1.10 negative            LAYOFF         EMPLOYEES  JOB REDUCTIONS         13     # 3     3    A  3.00 aneutral         POLLUTION           NUCLEAR    POWER PLANTS         13     # 4     4    A  1.40 positive CHEMICAL DISASTER             FUELS       POLLUTION         14     # 5     5    A -0.20 positive     PRESS RELEASE STOCK OPTION PLAN       EMPLOYEES         13     # 6     6    A  0.60 negative            PEOPLE        EXECUTIVES           FRAUD         14     # 7     7    A  0.10 aneutral         EMISSIONS               CO2  CLIMATE CHANGE         12     # 8     1    B -0.21 positive            ENERGY             SOLAR  SUSTAINABILITY         11     # 9     2    B -1.20 negative  WASTE MANAGEMENT         POLLUTION HAZARDOUS WASTE         13     # 10    3    B  0.90 negative         EMPLOYEES        EXECUTIVES       BONUS PAY         13     # 11    4    B  0.30 positive         SUBJECT11         SUBJECT21       SUBJECT31         14     # 12    5    B -0.10 aneutral         SUBJECT12         SUBJECT22       SUBJECT32         13     # 13    6    B  0.30 aneutral         SUBJECT13         SUBJECT23       SUBJECT33         14     # 14    7    B -0.12 aneutral         SUBJECT14         SUBJECT24       SUBJECT34         12      This regression model would look like this:               mymodel1 <- lm(RET ~ CLASS + CONTROLVAR, data=df)      aneutral (neutral) will be the reference category. I would also like to see the effect of certain subjects of the article on the abnormal return. How can I do that? Let's say I want to include the subjects LITIGATION, POLLUTION and LAYOFF. I'd like to see how the effect of positive, negative and neutral news change, if the article is about POLLUTION for example. If I make "dummy columns" for the three subjects of interest, my data.frame looks like this:               DATE <- c("1","2","3","4","5","6","7","1","2","3","4","5","6","7")     COMP <- c("A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B")     RET <- c(-2.0,1.1,3,1.4,-0.2, 0.6, 0.1, -0.21, -1.2, 0.9, 0.3, -0.1,0.3,-0.12)     CLASS <- c("positive", "negative", "aneutral", "positive", "positive", "negative", "aneutral", "positive", "negative", "negative", "positive", "aneutral", "aneutral", "aneutral")     SUBJECT.1 <- c("LITIGATION","LAYOFF","POLLUTION","CHEMICAL DISASTER","PRESS RELEASE","PEOPLE","EMISSIONS","ENERGY","WASTE MANAGEMENT","EMPLOYEES","SUBJECT11","SUBJECT12","SUBJECT13","SUBJECT14")     SUBJECT.2 <- c("POLLUTION","EMPLOYEES","NUCLEAR","FUELS","STOCK OPTION PLAN","EXECUTIVES","CO2","SOLAR","POLLUTION","EXECUTIVES","SUBJECT21","SUBJECT22","SUBJECT23","SUBJECT24")     SUBJECT.3 <- c("ENVIRONMENT","JOB REDUCTIONS","POWER PLANTS","POLLUTION","EMPLOYEES","FRAUD","CLIMATE CHANGE","SUSTAINABILITY","HAZARDOUS WASTE","BONUS PAY","SUBJECT31","SUBJECT32","SUBJECT33","SUBJECT34")     LITIGATION <- c(1,0,0,0,0,0,0,0,0,0,0,0,0,0)     POLLUTION  <- c(1,0,1,1,0,0,0,0,1,0,0,0,0,0)     LAYOFF     <- c(0,1,0,0,0,0,0,0,0,0,0,0,0,0)     CONTROLVAR <- c("11","13","13","14","13","14","12","11","13","13","14","13","14","12")          df2 <- data.frame(DATE, COMP, RET, CLASS, SUBJECT.1, SUBJECT.2, SUBJECT.3, LITIGATION, POLLUTION, LAYOFF, CONTROLVAR, stringsAsFactors=F)          df2           #    DATE COMP   RET    CLASS         SUBJECT.1         SUBJECT.2       SUBJECT.3 LITIGATION POLLUTION LAYOFF CONTROLVAR     # 1     1    A -2.00 positive        LITIGATION         POLLUTION     ENVIRONMENT          1         1      0         11     # 2     2    A  1.10 negative            LAYOFF         EMPLOYEES  JOB REDUCTIONS          0         0      1         13     # 3     3    A  3.00 aneutral         POLLUTION           NUCLEAR    POWER PLANTS          0         1      0         13     # 4     4    A  1.40 positive CHEMICAL DISASTER             FUELS       POLLUTION          0         1      0         14     # 5     5    A -0.20 positive     PRESS RELEASE STOCK OPTION PLAN       EMPLOYEES          0         0      0         13     # 6     6    A  0.60 negative            PEOPLE        EXECUTIVES           FRAUD          0         0      0         14     # 7     7    A  0.10 aneutral         EMISSIONS               CO2  CLIMATE CHANGE          0         0      0         12     # 8     1    B -0.21 positive            ENERGY             SOLAR  SUSTAINABILITY          0         0      0         11     # 9     2    B -1.20 negative  WASTE MANAGEMENT         POLLUTION HAZARDOUS WASTE          0         1      0         13     # 10    3    B  0.90 negative         EMPLOYEES        EXECUTIVES       BONUS PAY          0         0      0         13     # 11    4    B  0.30 positive         SUBJECT11         SUBJECT21       SUBJECT31          0         0      0         14     # 12    5    B -0.10 aneutral         SUBJECT12         SUBJECT22       SUBJECT32          0         0      0         13     # 13    6    B  0.30 aneutral         SUBJECT13         SUBJECT23       SUBJECT33          0         0      0         14     # 14    7    B -0.12 aneutral         SUBJECT14         SUBJECT24       SUBJECT34          0         0      0         12      The first problem is, that the dummy variables partially overlap. My model would look something like this.               mymodel2 <- lm(RET ~ CLASS + LITIGATION + POLLUTION + LAYOFF                      + LITIGATION * CLASS + POLLUTION * CLASS + LAYOFF * CLASS   # Interaction Variables                    + CONTROLVAR,     # Control Variables                    data=df2)      I'm quite sure this model is wrong, but I don't know exactly what is wrong and how to implement a model that works for this task. Can anyone help me with this problem? Thank You!


title: Multiple linear regression for hypothesis testing

text: I am familiar with using multiple linear regressions to create models of various variables. However, I was curious if regression tests are ever used to do any sort of basic hypothesis testing. If so, what would those scenarios/hypotheses look like?


title: Is Bayesian nonlinear regression using conjugate priors possible?

text: I have data for many time series. Each time series can be modeled by a particular nonlinear function of several parameters and, though I want point estimates of the model parameters for each of the many time series, I am more interested in the uncertainty associated with the parameter estimates (i.e., in the strength of the data). I have good prior information that I'd ideally like to include. I have tried MCMC, but am not confident that reasonable quality posterior samples can be obtained for the many time series (I can tweak the parameters of the MCMC algorithm for specific time series, but I'm not convinced my settings generalize well for other time series); additionally, CPU time is an issue. Rather than use MCMC, I'm wondering about using conjugate priors (a multivariate normal is a good model of my log-transformed parameters). However, I'm having trouble finding examples of Bayesian nonlinear regression with conjugate priors. Is it possible? Your assistance is much appreciated.


title: Trend analysis in temperature time series

text: Does anyone know how to determine the significance of the slope from linear model fitted using generalised least squares (GLS)? I am fitting a linear model to temperature time series with the aim of assessing for trends using gls function in nlme package, but I could not figure out the p-value (s) for the slope in model outputs. I chose GLS because of serial correlation.


title: Regarding using bigram (N-gram) model to build feature vector for text document

text: A traditional approach of feature construction for text mining is bag-of-words approach, and can be enhanced using tf-idf for setting up the feature vector characterizing a given text document. At present, I am trying to using bi-gram language model or (N-gram) for building feature vector, but do not quite know how to do that? Can we just follow the approach of bag-of-words, i.e., computing the frequency count in terms of bi-gram instead of words, and enhancing it using tf-idf weighting scheme. Thanks.


title: How to interpret a logistic regression model with all negative coefficient?

text: I have 4 predictors, and 1 binary response. I fitted a logistic regression model. A strange thing is that all the coefficient of the model are negative. Is that possible? Probably I did something wrong. My interpretation is that the odds ratio of either variable is less than 1. That is, neither variable actually do any good to the response. Even the intercept is negative. Please share your intelligence.


title: How to get asymptotic covariance matrix when observed information matrix is singular

text: I'm fitting different models by Maximum Likelihood. To do this I'm using a stochastic version of Newton-Raphson algorithm, where both the gradient and the Hessian of the likelihood are estimated at each step. In most cases I reach convergence, but then I often encounter the following problem: the estimate Hessian $\hat H$ at convergence is negative definite. This is a problem because I can't invert it to get an estimate of the observed Fisher information. This happens because one or more parameters are weakly identified. On the other hand other parameters seem to be well identified and their corresponding entries in the estimated Hessian seem to make sense. Identifiability is difficult to assess beforehand for the dynamical models I'm trying to fit. What I would like is an approach that (starting from the $\hat H$) points out what set of parameters are not identifiable and that gives variances for the remaining parameters. I tried to do this by tilting the smallest eigenvalue of $\hat H$ in order to get a better conditioning number. The results seem arbitrary: depending on the conditioning number I want to achieve I get different variances for the parameters. EDIT: This is a typical example of an Hessian I can end up with:               H <- matrix(c( 67.23586, 10.477815, 138.696877,                    10.47782, -3.238982,   2.592774,                   138.69688,  2.592774, 473.161347 ), 3, 3, byrow = TRUE)      You see that I have a negative entry in the main diagonal: the second parameter is weakly identifiable, while the other well identified. The other situation I often encounter is that of sub-group of parameters that are redundant - highly correlated.


title: standard error of known population values

text: i tried to get a clearer understanding of the standard error by constructing a hypothetical population consisting of the following values: 1, 3, 5, 7. i calculated the sample means of all samples with a size = 2. as expected, the mean of all sample means was equal to the population mean (4). however, the standard deviation of the sample means (1.633)--also known as the standard error of the mean was not equal to the population standard deviation (2.58) divided by the square root of the sample size (square root of 2) which was 1.825. but, i discovered that the variance of the sample variances was 5.33 and the standard deviation of the samples variances was 2.31. dividing 2.31 by the square root of the sample size (square root of 2) yielded a value of 1.633 which is identical to the standard deviation of the sample means. i am bit confused because the standard error according to most authors is the population standard deviation divided by the square root of the sample size. however, in this case the standard error is equal to the standard deviation of the sample variances divided by the square root of the sample size. i would appreciate if you could help me on this matter. thank you.


title: Box-Cox transform units and scaling

text: Why is the Box-Cox transform defined the way it is, in particular with regards to the geometric mean? In particular, if you're paying attention to units, wouldn't you need to use $\bar{y}\log(y/\bar{y})$ for $\lambda=0$, and $\bar{y} \cdot ((y/\bar{y})^\lambda - 1)/\lambda$ for $\lambda \ne 0$, with $\bar{y} \equiv \mathrm{GM}(y)$ to avoid either taking the log of or subtracting 1 from a value with dimensions? What would be the advantage of the "standard defintion" over this variant?


title: Question about logistic regression

text: I want to run a binary logistic regression to model the presence or absence of conflict (dependent variable) from a set of independent variables over a 10 year period (1997-2006), with each year having 107 observations. My independents are:   * land degradation (categorical for 2 types of degradation);    * population increase (0- no; 1-yes);    * livelihood type (0 - type one; 1 - type two);    * population density (three levels of density);    * NDVI continuous (max. veg productivity);    * NDVI$_{t-1}$ (decline in veg from the previous year - 0 - no; 1 -yes) and   * and NDVI$_{t-2}$ (decline in veg from two years past - 0- no; 1- yes).  I am fairly new to it all - this is a project my lecturer has given me - and so I would be grateful of some advice or guidance. I have tested for multicolliniarity already. Essentially my data is split up into 107 units of observation (spatial regions) covering 10 years (1070 in total) and for every unit of observation it gives be a 'snapshot' value of conditions of the independent variables at that time within that unit (region). I want to know how to set up my logistic regression (or table) to recognize the 107 values of each year separately so that the temporal NDVI changes between different unit years can be assessed?


title: Probability distribution of windowed cross-correlation

text: This question is in the context of time-delay estimation. Say I have a stationary Gaussian stochastic process $g$, and I know its autocorrelation function $R_g(\tau)$. To do time-delay estimation, I'm computing a windowed cross correlation between $g$ and a delayed version of it. In other words, $$ g_1 = g(x-D) \\\ \phi(\tau) = \int_{-T/2}^{T/2} g(x) g_1(x + \tau) $$ and I'm going to determine the delay by finding the maximum of $\phi$. My question is, is it possible to get an expression for the probability distribution of $\phi$?


title: Forecast with arima model

text: Starting with arima models in R, I cannot make a forecast with my favourite model. For example, the commands               predict(arima(data_ts,order=c(1,1,2),xreg=cbind(t),seasonal=list(order=c(0,1,1),period=12)))       and               forecast(arima(data_ts,order=c(1,1,2),xreg=cbind(t),seasonal=list(order=c(0,1,1),period=12)))      don't work. Can you explain me why?


title: Calculating Standard deviation of percentages?

text: I have the following data X 1 2 3 4 5… Y 10 12 13 14 15… X/Y 10% 16% 23% etc. How do I find the standard deviation of percentage (last line)? Can I treat the ratio as a normal distribution and apply regular SD formula?


title: Convergence failure with step-by-step parameter estimates in SAS proc logistic using the OFFSET option

text: I'm trying to fit a logistic model of the following form: $$ y \sim \frac{1}{1 + \exp(-\beta_X)} $$ where $$ \beta_X = a_0 + a_1x_1 + a_2x_2 + a_3x_3 + a_4x_4 + a_5x_5 $$ In my case, there is often good reason to first fit to a simpler model of $x_1$ and $x_2$; After that, fix $a_0, a_1$ and $a_2$, and then fit the rest. With PROC LOGISTIC, one can achieve that by doing the following:               DATA myData;        SET myData;       restrict = a0 + a1*x1 + a2*x2;       RUN;          PROC LOGISTIC DESCENDING DATA=myData;       MODEL y = x3 x4 x5 / LACKFIT NOINT OFFSET=restrict;       RUN;      The trick is to use the OFFSET option. That's all fine. The problem is, one often encounters convergence failures. Below is a typical log from SAS.               WARNING: Ridging has failed to improve the loglikelihood. You may want to use a different              ridging technique (RIDGING= option), or switch to using linesearch to reduce the              step size (RIDGING=NONE), or specify a new set of initial estimates (INEST=              option).           WARNING: The LOGISTIC procedure continues in spite of the above warning. Results shown are              based on the last maximum likelihood iteration. Validity of the model fit is              questionable      Does anyone have experience with such problems? I'm looking for (1) an understanding of why the issue occurs (2) suggestions on how to get around it. One question comes to mind, does constraining $a_0, a_1, a_2$ change the globally concave property of the maximum likelihood function of the logistic model with respect to the remaining list of parameters: $a_3, a_4, a_5$? If it doesn't, why is there a convergence issue at all?


title: How to calculate the missing values of this sequence effectively?

text: We have a sequence $\lbrace b \rbrace$ which is defined in a relation $b_i = b_{i-1} * 0.9 + t_i * 0.1$. If we have observations of the full sequence of $\lbrace b \rbrace$, it would be quite simple to calculate the values of $\lbrace t \rbrace$. However, the case is that we only have around $50\%$ to $80\%$ observations in $\lbrace b \rbrace$ and we still want to calculate the values of $\lbrace t \rbrace$ as accurately as possible. Are there any well- known mathematical methods to solve this problem? Thanks!


title: Distribution of inverse Wishart to a power?

text: In a related question, I had asked about the norm induced by an inverse Wishart matrix. I am interested in generalizing that result somewhat. Let $A\sim\mathcal{W}_p\left(I,n\right)$, a Wishart matrix with scale matrix $I$, the identity, and $n$ degrees of freedom. Let $l$ be some fixed $p$-dimensional vector. Consider $$h = \frac{l^{\top}l}{l^{\top} \left(A^{-1}\right)^{m} l},$$ where $m$ is some integer. When $m=1$, $h$ is evidently distributed as a Chi-square random variable (with $n-p+1$ degrees of freedom, I believe). Are there known results on the distribution of $h$ when $m=2$? Other values?


title: How to detect regression outliers when many residuals are zero

text: I have a situation that when I fit a robust regression line (least trimmed squares) to a set of data a lot of the residuals are in fact zero. This occurs mainly in the situation where the slope is zero and the y values are integers. When the line is fitted it runs right through the majority of the values. I am happy with this fit as a traditional least squares line is incorrect due to a few errors in the data. However now I've fitted the line I want to analyze the residuals and detect the outliers (possibly automatically). I intended to compute a 'score' perhaps based on Tukey's technique of using the upper and lower quartile plus 1.5 times the inter-quartile range. However this approach won't work if many of the residuals are zero because the IQR is also zero. What should I do? Just base the score on something like the number of standard deviations from the mean? Thanks.


title: How to obtain mean and joint distribution in factor analysis?

text: I have a factor analysis model defined by: $x = m + Wz + e$ where $x$ is a p-dimensional visible variable, $m$ is a constant vector, and $z$ is a $n$-dimensional Gaussian latent variable with $z$ ~ $N(0, I)$, $W$ is a $p\times m$ matrix and $e$ is a $p$-dimensional with $e$ ~ $N(0, \Psi$) - it is the factor loadings matrix. $Z$ and $e$ are independent. $p(x)$ is defined by the model is Gaussian, but **how do I find its mean?** **Also, what is the explicit joint distribution of $p(z, x)$?**


title: How do I fit a nonlinear mixed effects model for repeated measures data using nlmer()?

text: I am trying to analyze repeated measures data and am struggling to make it work in `R`. My data is essentially the following, I have two treatment groups. Every subject in each group is tested everyday and given a score (the percentage correct on a test). The data is in the long format:               Time Percent Subject   Group        1       0    GK11 Ethanol        2       0    GK11 Ethanol        3       0    GK11 Ethanol        4       0    GK11 Ethanol        5       0    GK11 Ethanol        6       0    GK11 Ethanol      The data resembles a logistic curve, subjects do very poorly for a few days followed by rapid improvement, followed by a plateau. I'd like to know if the treatment has an effect on the test performance curve. My thought was to use `nlmer()` in the `lme4` package in `R`. I can fit lines for each group using the following:               print(nm1 <- nlmer(Percent ~ SSlogis(Time,Asym, xmid, scal) ~ Asym | Subject,     salinedata, start = c(Asym =.60,  xmid = 23, scal = 5)), corr = FALSE)      I can than compare groups by looking at the estimates for the different parameters and standard deviations of the estimated lines but I'm not sure this is the proper way to do it. Any help would be greatly appreciated.


title: The tanh activation function in backpropagation

text: In the backpropagation algorithm when the output activation function is `tanh` and the number of classes is 2 (binary problem), the value obtained at the output layer is in the range between -1 to 1. The cross-entropy error function has `log` that is applied on the predicted values. Therefore, if one of the output values is a negative number, an invalid operation, namely, `log (non- positive number)`, occurs, rendering the cross-entropy function invalid. This boils down to the following questions:   * Is it disallowed to set the output activation as `tanh`?    * Should the output activation always be the softmax even for a binary class   problem?


title: Conditional on Gaussian, need clarification

text: I'm reading Andrew Ng's notes on machine learning, and on page 12 of this document, he makes a step in his proof that I'm trying to decipher: Let $\textbf{x} = \left( 1 , x_1 , x_2 , \cdots , x_n \right)^T$, a vector of variables, and $\theta = \left( \theta_0 , \theta_1 , \theta_2 , \cdots , \theta_n \right)^T$, a vector of linear coefficients of those variables. Let's define $y$ as $$y = \theta ^T \textbf{x} + \epsilon$$ where $\epsilon \sim \mathcal{N}(0,\sigma^2)$, that is $$p(\epsilon) = \frac{1}{\sqrt{2 \pi}}\exp \left( -\frac{\epsilon^2}{2\sigma^2}\right).$$ Next line says the following about conditional probability of $y$ given $\textbf{x}$ and coefficients $\theta$, which are treated as deterministic: $$p(y|x;\theta) = \frac{1}{\sqrt{2 \pi}}\exp\left( -\frac{(y- \theta ^T \textbf{x})^2}{ 2\sigma^2}\right)$$ Can someone help me see how we get this conditional distribution?


title: Where did this risk exposure 'estimation-formula' come from?

text: I was reading a book and the authors metioned that risk exposure can be estimated _scientifically_ using this forumula: $risk(\$) = \frac{(a + 4m + b)}{6}$ and standard deviation $\sigma = \frac{b-a}{6}$ Where:   a = minimum dollar exposure   m = most likely dollar exposure   b = maximum dollar exposure   Where `a, m, b` are solicited from past data or expert judgement. All this is fine, but I've never seen this formula - where did this come from? Is it founded in statistics or any scientific reference? Anyone know anything about this? Pointers/suggestions/clarifications would be greatly appreciated. I'm curious as to how does the above formula provides a risk exposure metric. Why '6' in the denomiator and why the `a + 4m + b` in the numerator? Similarly for $\sigma$? Any ideas? **EDIT** : This formula seems to be that of weighted averaging but am not sure about it's inception or the reason of giving 'm' a weight of 4 and how that influences the calculation of $\sigma$? **UPDATE** :Book - Making the Software Business Case. However just Googling for the numerator threw up a few results - it seems to be known as the PERT formula and is used for estimation of 'time' - but my question about the 'formula' still holds irrespective of it being used for risk... **DERIVATION OF FORMULA** \- The reason for the formula to be what it is seems to be beautifully explained in a JSTOR article for those interested!


title: Looking for an efficient algorithm to detect Tomek links

text: I'm looking for an efficient algorithm to detect Tomek links. I'm wondering if anyone knows where to find it. This is the definition of Tomek links: Suppose $\\{ E_1,\ldots,E_n\\} \subset R^k$ is a dataset, with each $E_i$ having exactly one of two labels $+$ or $-$. A pair $(E_i,E_j)$ is called a Tomek link if $E_i$ and $E_j$ have different labels, and there is not an $E_l$ such that $d(E_i,E_l) < d(E_i,E_j)$ or $d(E_j,E_l) < d(E_i,E_j)$, where $d(x,y)$ is the distance between $x$ and $y$. Thanks!


title: Belief propagation on MRF with complex cliques

text: Is there a belief propagation algorithm for exact inference on a MRF with complex clique structures (i.e. ones involving more than 2 neighbours)? For MRF's with cliques that only involve pairwise interaction, you could just search out far enough and cluster to form an acyclic graph and run the usual BP. With more complex cliques, this seems impossible to me as clustering might involve cutting through a clique with multiple members on either side. Is there a workaround for this? Perhaps some clever conditioning arguments?


title: Feature selection for the text mining?

text: Before performing the task of text mining, we need to select the features for characterizing each given document. Are there any systematic guidance on choosing the document features? How does the length of the document affect the feature selection process for the documents?


title: How to fill in missing values to a data set?

text: I have a picture from a very old archive. It is a list of tuples (var, percentage). The problem is that I have only the 13 top values, ordered on percentage. In python I decided to write it as a dictionary. myolddictionary = {'var1': 0.64, 'var2': 0.63, 'var3': 0.40, 'var4': 0.38, 'var5': 0.36, 'var6': 0.34, 'var7': 0.34, 'var8': 0.34, 'var9': 0.33, 'var10': 0.32, 'var11': 0.32, 'var12': 0.30, 'var13': 0.29} Now, my intuition tells me that if I had the entire dataset, because of the very small percentages of the top values, it would look like a uniform distribution, or a normal distribution with a very high variance, because the number of vars can be more than 100. My initial thought was to just take sum of the percentages I have, take the difference from 100, and then divide it into that many new vars so that each one gets a percentage smaller than 0.29% (otherwise, they would appear in the archive, right?) and add them to the dictionary. But then I thought that this might be a really stupid way. What I actually want is to draw from a distribution where: a)draws are not equally probable b)the sum of all draws is what remains from 100 after removing the percentages in my initial list (95.01 in this case) c)all draws are smaller or equal to 0.29. Something like function superamazingdistribution(number_of_extra_draws = 0, sum_of_draws, range) where: 1)number_of_extra_draws is the number of extra draws I want to do, added to the minimum required so can reach the sum_of_draws using the minimum number of draws * range 2)sum_of_draws, the number I want to reach, in this particular case 95.01 3)a list [num1, num2], in my case [0, 0.29] and it returns a dictionary (one of many possible I guess) like {var1: percentage1, var2: percentage2,...,varN: percentageN}, which I can loop through and add to my own dictionary. Sorry if this is a stupid question, sadly I never took a statistics course in my life.


title: Estimating Poisson waiting times: using the last, censored sample?

text: Assume I have a process that I assume to be Poisson, meaning that the waiting time between events is exponentially distributed. As a concrete example, assume I want to estimate the waiting time to a car accident at an intersection. These accidents happen infrequently, so the data I have must be used efficiently. It takes $X_1$ days for the first accident to occur, $X_2$ days after the first accident for the second, etc. It's been $X_C$ days (and counting) since the last accident. What's the best (unbiased, minimum variance) way to estimate the waiting time in this situation, where it's been some time since the last Poisson event occurred and I want to take advantage of my knowledge of $X_C$? Edit: I am interested in the exponential rate parameter. The interesting part of the question is including the censored value $X_C$ in the estimate. The terminology for this somehow escaped me when I wrote the original question.


title: How to do post-hoc tests with logistic regression?

text: I looked at whether participants recycled after being subjected to a given message. There were three groups: a Control group, group A, and group B. Since my dependent measure was dichotomous (recycling vs. using the waste bin), I performed a logistic regression in SPSS. However, SPSS only let me compare individual groups to the Control group. I found that group A recycled more than controls, but group B did not. What I am curious about now is whether group A recycled significantly more than group B. Is it possible to do such a post-hoc test? If so, how?


title: What do you think is the best goodness of fit test?

text: I am looking at fitting distributions to data (with a particular focus on the tail) and am leaning towards Anderson-Darling tests rather than Kolmogorov- Smirnov. What do you think are the relative merits of these or other tests for fit (e.g. Cramer-von Mises)?


title: Data Mining-- How to Tell Whether the Pattern Extracted is Meaningful?

text: I am sure that everyone who's trying to find patterns in historical stock market data or betting history would like to know about this. Given a huge sets of data, and thousands of random variables that may or may not affect it, it makes sense to ask any patterns that you extract out from the data are indeed true patterns, not statistical fluke. A lot of patterns are only valid when they are tested in the samples. And even those that are patterns that are valid out of samples may cease to become valid when you apply it in the real world. I understand that it is not possible to completely 100% make sure a pattern is valid all the time, but besides in and out of samples tests, are their any tests that could establish the validness of a pattern?


title: Open source tools for visualizing multi-dimensional data?

text: Besides gnuplot and ggobi, what open source tools are people using for visualizing multi-dimensional data? Gnuplot is more or less a basic plotting package. Ggobi can do a number of nifty things, such as:   * animate data along a dimension or among discrete collections   * animate linear combinations varying the coefficients   * compute principal components and other transformations   * visualize and rotate 3 dimensional data clusters   * use colors to represent a different dimension What other useful approaches are based in open source and thus freely reusable or customizable? **Please provide a brief description of the package's abilities in the answer.**


title: What are the distributions on the positive k-dimensional quadrant with parametrizable covariance matrix?

text: Following zzk's question on his problem with negative simulations, I am wondering what are the parametrized families of distributions on the positive k-dimensional quadrant, $\mathbb{R}_+^k$ for which the covariance matrix $\Sigma$ can be set. As discussed with zzk, starting from a distribution on $\mathbb{R}_+^k$ and applying the linear transform $X \longrightarrow\Sigma^{1/2} (X-\mu) + \mu$ does not work.


title: Dichotomous variables analysis

text: I am analyzing data from different dichotomous variables. When analyzing association between these variables, what is the most recommended method?


title: Gradient boosting in R uses only a single variable

text: I am trying to build a boosting model using the package gbm in R. I have the following code:               gb = gbm(aaa_target ~ .,          data=myDdata,          n.trees=100,          verbose=TRUE)      and when I have trained the model, I can get a summary like this:               summary(gb)      The issue I am having, is that only a single variable (out of around 30) is selected and is given 100% predictive power. I know for a fact that many of the variables carry information (although the selected one is the most significant one), and using the randomForest package gives me a model which assigns significance to many of the variables. Does anybody have a clue to why this might be the case?


title: Converting 2nd order Markov chain to the 1st order equivalent

text: Given a 2nd order Markov chain where each state takes values in the set $\mathcal{X}=\\{A,C,G,T\\}$, such that all transition probabilities $p(x_t|x_{t-1},x_{t-2})$ are larger than zero, How to convert it to the equivalent 1st order one with all the transition probabilities defined?


title: Regression where a subset of observations are missing data on an independent variable

text: Consider the regression equations below: \begin{align} Y_i &= \beta_0 + \beta_1 X_{i1} + \varepsilon_i \\\ Y_j &= \beta_0 + \beta_1 X_{j1} + \beta_2 X_{j2} + \varepsilon_j \end{align} where $Y_i,\ X_{i1},\ \varepsilon_i,\ Y_j,\ X_{j1},\ \& \ X_{j2},\ \varepsilon_j$ are vectors, and $_i$ and $_j$ index distinct sets of observations. The $_i$ respondents did not meet a qualification criterion and hence were not asked the question that corresponds to $X_2$. The dependent variable and the first independent variable is the same in both regression equations but the second regression equation has an independent variable that is not present in the first. Obviously, I can estimate the two regressions separately but that will not be efficient. Therefore, I was considering re-writing the first one as: $$ Y_i = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \varepsilon_j $$ where $X_{i2}$ is a vector of $0$s. Then I can estimate the parameter estimates by using OLS with the equation below: $\left[ \begin{array}{ccc} Y_i\\\ Y_j\end{array} \right] = \left[ \begin{array}{ccc} {\bf 1} & X_{i1} & X_{i2} \\\ {\bf 1} & X_{j1} & X_{j2}\end{array} \right] \left[ \begin{array}{ccc} \beta_0\\\ \beta_1\\\ \beta_2\end{array} \right] + \left[ \begin{array}{ccc} \epsilon_i\\\ \epsilon_j\end{array} \right]$ In the above equation, ${\bf 1}$ stands for a vector of $1$s of the appropriate dimension. Is the above a standard approach to obtaining efficient estimates? Is there a name to this way of estimation?


title: Weighted Random Selection with bounded replacement

text: The problem: Let there be an array (noted as X) where each item has a member, and a label assigned to it. For example: $$ X= \\{a_0, b_0, b_1, c_0, c_1, a_1\\} $$ Here the letter denotes a label, and the index is just another instance from the given label. Now the goal is to sample n items from X into Y, where: $$ |Y|=n < m = |X| $$ and each label has a weight assigned to it, so that in the long run the weights will define the distribution of the elements in Y. Elements from within a label are different, but whichever of them may be selected. The multiplicity of the labels are random too (it may vary from 1 to 10). The length of X will be usually around 100, from which around 10 items need to be selected. So, let's say if the weights are $$ W = \\{a=5, b=45, c=50\\} $$ for 10,000 runs if we add up the elements from Y it will be true that 5% of it will be from a, 45% from b, and 50% from c. The indexes just denote that you have multiple item instances from a label. There is no frequency requirement for within label (so simple random sampling can/should be used). The sampling is without replacement (considering the full array (X), however it is with bounded replacement if you consider just the labels (that is a,b,c labels). At the end we do not need exact match with W, expected value is okay. Successive runs may result in slightly different output, that too is okay (although the expected value of the frequency of labels, should approximate W). Now, I've already tried to do this in python via the roulette wheel selection, where if a label has no more items I redo the roulette, however this sometimes goes really off.   1. What's a good algorithm for this?   2. What's wrong with my code? (why the big difference sometimes - here's the python code http://pastebin.com/SFhV74Z8)


title: Finding probability density function with unknown values

text: I am not sure if this is the place to ask, have tried to read up on probability density function (PDF) in order to answer this question but to no avail. How do I go about starting on this? How do I generate a function based on the given $\lambda$ and that small $t$ symbol. What is $X \sim Ga(n,\lambda)$? An explanation before the answer would be nice. Appreciate all the help needed. Regards ![Stats Question](http://i.stack.imgur.com/OzD7G.png)


title: survival analysis without enough data

text: I have computed and plotted the survival function for a subscription-based service and the following is the result. ![enter image description here](http://i.stack.imgur.com/qm7dm.png) The problem is that there does not seem to be enough data to get a full curve. This is because most of the oldest accounts are still active. So my question is would it still be useful to compare survival curves for different segments of accounts given that there is clearly not enough data for a full curve.


title: Correlogram in R like in Stata?

text: In STATA I can create a "Correlogram" to find the appropriate lag order in case of time series. E.g. ![corrgram](http://i.stack.imgur.com/8rt9V.png) I know I can use the acf or Acf of the forecast package to calculate the ACF and PACF and to plot it. But how can I get the significance values? I mean the values in the column Prob>Q? Is this implemented in any package/command?


title: Is this an ARMA(1,1) or something else?

text: I have a company, for which I calculated the logarithmic loss (-logarithmic return). Now I want to fit a mean equation to the returns, so I have to think about fitting an ARMA(p,q). I lookt at the acf and pacf and get the following picture: ![arma](http://i.stack.imgur.com/nvXop.png) This seems not to be a simple AR or MA? What is this, what should I use? Data can be found here. It would be nice, if you explain me what you do in order to know what process is appropriate.


title: Omit 0 lag order in ACF plot

text: How can I omit the zero lag order in an acf plot? See this picture: ![enter image description here](http://i.stack.imgur.com/RP3u1.png) generated by               dummy<-c(14,0.004,0.2,1,0.002,-3,-0.042,1.2,-1,1.3,2.1,4,3001,-2,0.3,2,3)     acf(dummy)      The high peak (which is logically 1) is destroying the plot, since the scaling is too big. I would like to omit the high peak at lag order 1, so that the scaling can be reduced to -0.2 up to 0.2 for example, how can I do this?


title: Finding Degree of Freedom from 2 Logit Models

text: Hi I have a question on finding DF in logit models. Say I have a table (all fictional numbers) ![8 cities smoking vs heart disease](http://i.stack.imgur.com/ZQxvw.png) and say I am testing goodness of fit (deviance) between the fitted model and saturated model. 1) main effect model : $$Logit[\pi (x,y))] = \alpha + \beta_{x}^{smoking} + \beta_{y}^{city}$$ 2) Saturated model: $$ Logit[\pi (x,y))] = \alpha + \beta_{x}^{smoking} + \beta_{y}^{city} + \beta_{xy}^{smoking,city}$$ How do I calculate the DF from these 2 models? (the answer is 7 from stats software, but I don't understand how to get the answer manually). DF=degree of freedom = difference in number of parameters from the 2 models. thanks


title: Linear SVM decision boundary after a linear transformation of data

text: Let $w$ be the decision boundary of a linear SVM trained on the dataset $D=\\{(x_i, y_i)_{i=1}^N\\}$. Suppose we apply a linear transformation A to examples $x_i$s and obtain a new dataset $D'=\\{(z_i, y_i)_{i=1}^N\\}$ where $z_i = Ax_i$. And let $w_2$ denote the decision boundary of the linear SVM trained on $D'$. Can we say anything about the relationship between $w$ and $w_2$? Depending on the properties of $A$, would it be possible to cheaply obtain $w_2$ from $w$?


title: How to visualize a 3-D density function?

text: What is the best way graphically to visualize a 3-D density function? As in I would like to visualize $$z=f_{X,Y}(x,y)$$? Not necessary but `R` code for this would be great.


title: Reporting regression equations for non-significant results

text: Do you report the regression equation for non-significant multiple regressions if you have reported them for significant regressions in the same work? My paper has six multiple regressions in it. Three are significant and I have reported the regression equation for each of these (e.g., $y=b_1+b_2x+a$); I can't find out if I need to report the equation for non-significant results too. Would it be best to just report them all or are there rules about it?


title: Looking for a good and complete probability and statistics book

text: I never had the opportunity to visit a stats course from a math faculty. I am looking for a probability theory and statistics book that is complete and self-sufficient. By complete I mean that it contains all the proofs and not just states results. By self-sufficient I mean that I am not required to read another book to be able to understand the book. Of course it can require college level (math student) calculus and linear algebra. I have looked at multiple books and I didn't like any of them.   * DeGroot & Schervish (2011) _Probability and Statistics (4th Edition)_ Pearson This is not complete enough. It just states a lot of stuff without the derivation. Besides that I like it.   * Wasserman (2004) _All of Statistics: A Concise Course in Statistical Inference_ Springer. Didn't like it at all. Almost no explanations. "Weighing the Odds" from David Williams is more formal than DeGroot and seems to be complete and self-sufficient. However, I find the style strange. He also invents new terms that only he seems to use. All the stuff that is explained in DeGroot too is explained better there. If you know a great book in German that's also fine as I am German.


title: Finding the distribution of a "row-transition"-based quantity computed over row-permutations of a matrix

text: Say we have N observations of triplets of binary variables V1,V2,V3. (i.e. one such observation might be: [V1 V2 V3] = [1 0 1] ). Now suppose we wish to arrange these observations in a random sequence/ordering (observation1; observation2; ....; observationN). Any such ordering can be represented as an Nx3 matrix, where the i-th row represents the i-th observation. For example, one ordering of these N rows would be: \begin{array}{ccc} 1 & 0 & 1 \\\ 0 & 0 & 1 \\\ ... \\\ 0 & 1 & 0 \end{array} Since there are 3 variables, and all variables are binary, there are 8 possible distinct "types" of observations (i.e. row vectors such as [0 1 1], [1 0 1], etc). For each of these types, suppose _we know how many rows (out of the N total rows) are of this type._ Now, suppose that for each possible row-permuted (Nx3) matrix (i.e. each possible ordering of the N observations) we compute a quantity H as follows:   * If, for any adjacent pair of rows, row i and row i+1, of the matrix, row i+1 has a different value for V1 than row i, but the same value of V2 as row i, then within this row-pair:      * For the row with V1=0, subtract the value of V3 in this row from the value of V3 in the row where V1=1 . For example, if we have a pair of adjacent rows: \begin{array}{ccc} 1 & 0 & 0 \\\ 0 & 0 & 1 \end{array} then we would compute (0 - 1) = -1 for this pair. We would find the same result if the rows were in the opposite order: \begin{array}{ccc} 0 & 0 & 1 \\\ 1 & 0 & 0 \end{array}   * Let H be the **sum of all such computed values** within a given row-permuted matrix, divided by the number of adjacent pairs of rows (with the property just described) that appear in the matrix. **My goal:** I am looking for a way to find the probability distribution (or at least the mean or median) of values of H (over all possible orderings of rows) analytically, or at least in some way that might yield a greater degree of 'understanding' than simply sampling from a computer simulation. Any help, insights, and advice on how one might begin approaching such a problem will be much appreciated! Perhaps, for example, there are other ways of describing and thinking about this problem that will lend themselves better to a (relatively) simple solution.


title: repeated measures ANOVA with covariate

text: I conducted a repeated measures ANOVA and added a continuous covariate. Now there is a significant main effect of that continuous covariate. If I report that main effect I need to say something about its direction... but how can I interpret its direction? The F statistic isn't very helpful here and I don't want to do a median split. Can SPSS report a B, beta or t statistic for the covariate?


title: Finding t distribution values for uniform and exponential functions

text: I know that a basic one-sample $t$ statistic can be calculated by $$t = \frac{x - \mu}{s / \sqrt{n}}$$ I also know that the following 4 functions $x\sim\mathcal{U}(0,1)$, $x\sim\mathcal{Exp}(2)$ (the exponential function where $\lambda=2$), $P(x=1/2) = 1$, and the set of $\\{P(x=025)=0.8, P(x=1.5)=0.2\\}$, all have an $\mathbb{E}(x)=0.5$. I've written a program to simulate the $t$-statistic for all of them, and have the values from $0, 0.25, 0.5,\ldots ,1.5$ with the values in 0.25 increments on the $t$-statistics for $n=10000$. I'm looking for a way to theoretically prove the following 2 equations $P(x=1/2) = 1$ and the set of $\\{P(x=025)=0.8, P(x=1.5)=0.2\\}$. The simulation values for the first all 1's from 0.5 to 1.5 (as would be intuitive here) and I want to express it theoretically. The values for the second are $0,0,0,0.8045,0.8045,\ldots ,0.8045$ for the set. I'm trying to theoretically reason out and mathematically express the proof and equations that express how these values are gotten (not simulated but calculated). For example, I'm trying to calculate the value for the set function and reason out why the simulated value is slightly above 0.8 at 0.8045 while the calculated value should be straight 0.8 from 0.5 to 1.5. I feel like I'm missing something that would be obvious to another set of eyes, but that I'm missing after even the 100th time I look at it over and over. **Edit:** What I'm mainly not understanding here is the set case. I'm trying to express why the simulated value is 0.8045 (just above 0.800) whereas the $P(x=0.25)=0.8$ which I would think would give a sim value of 0.800 but doesn't. I edited to make this clearer, but I'll re-read it a few more times to see if I can make it clearer still.


title: Sampling weight, logistic regression and chi-square analysis

text: I have some questions. I am doing logistic regression. but my sample is small(n=129). And in DV, the smallest number of cases is 27(1'). Actually, I hope to use 12 predictors. However, I know the rule of thumbs(1 to 10), too. Then, if I use sampling weight variable, can I use the 12 predictors? when I use weight variable, my sample changes large data(over 20,000). Additionally, when I am conducting the chi-square analysis through a same sample, should I use the weight variable which is sampling weight about the population. please share your opinion. thank you very much!!


title: restricting range of values for AR(1) process

text: We have a stationary AR(1) process defined as $A_{t}= d+ \rho A_{t-1}+ \epsilon_{t}$ with $d>0$; $-1<\rho<1$; $\epsilon_{t}$ is white noise and follows $N(0,\sigma^2)$. What conditions do we need to ensure that $A_{t}$ takes values in $(0,1)$ for all $t$?


title: Can effective sample size in MCMC simulation be greater than the actual sample size?

text: I used coda package's `effectiveSize()` to find the effective sample size of my MCMC simulation. My effective sample size is greater than the actual sample size, e.g. 9813.626 greater than 9501. I wonder if that makes sense? My understanding is that the effective sample size cannot be greater than the actual sample size, and will be smaller when there is more autocorrelation.


title: Dependent variable coding for logistic regression in R

text: Given the outcome variable in a dataframe is a factored/categorical variable, when regressing the dependent variable (DV) onto a set of independent variables (IVs), what is the model predicting? The probability that the DV is the first level of the factor? Or the second? A related question: I know that given a numerical column of $1$s and $0$s, a logistic regression would model the probability of the higher order variable (i.e., value=$1$), so I have been attempting to recode the factor "character" variable into a numerical one. I am coming from a SAS background, so I am entirely to used to `if var = "yes" then var_num = 1; else var_num=0;` That's clearly wrong. What's the most efficient way you have found to recode such variables?


title: How to handle large .csv file in R?

text: I have a large(>100,000) single column floating point time-series data. I want to find structural changes within the data with respect to time( in my case index). In-order to do that, I am using R strucchange package. But I'm getting error. I expect to generate the following kind of plot with my dataset ![enter image description here](http://i.stack.imgur.com/sU31W.png) Here is the code which I'm using:               data = read.csv("tst1.csv", header = FALSE)     library(strucchange)          ## fit, visualize and test OLS-based CUSUM test     ## with a Cramer-von Mises functional     ocus <- efp(data ~ 1, type = "OLS-CUSUM")     sctest(ocus, functional = "meanL2")          ## estimate breakpoints     bp <- breakpoints(data ~1)     summary(bp)          plot.ts(data)     lines(bp)     lines(fitted(bp, type = "mean"), col = 2)      * * * I'm getting following error messages:               Loading required package: zoo          Attaching package: ‘zoo’          The following objects are masked from ‘package:base’:          as.Date, as.Date.numeric          Loading required package: sandwich     Error in model.frame.default(formula, data = data) :      invalid type (list) for variable 'data'     Calls: efp -> model.frame -> model.frame.default     Execution halted      * * * My data tst1.csv looks like:               0.6     0.78     0.54     0.96     0.43     1.43     1.03     0.83     0.68     3.23     3.09     2.87     ...      * * * Any help is much more appreciated. Thanks


title: Compare regression slopes of repeated measures linear regression

text: In my design, I have two groups of subjects and every subject is tested in four different conditions. So, I have a within-subject factor ('span_num', which ranges from 0 to 3) and a between-subject factor (group, which can be 'Linear' or 'U-shape'). My goal is to show that the slope between the spans (from 0 to 3) is higher in the Linear group than in the U-shape group. The slopes look very different and are significantly different when I compare the regression models I get when I ignore that my within-factor is a within- factor and treat it as a between-factor. I compared the slopes like this (but I don't trust the comparison because I don't trust the SE values of the slopes because I am treating my within-factor like a between-factor):               linear_lm <- lm(RT ~ span_num, dat = data[data$group == "Linear",])     ushape_lm <- lm(RT ~ span_num, dat = data[data$group == "U-shape",])     linear_intercept <- summary(linear_lm)$coefficients[[1]]; linear_ise <- summary(linear_lm)$coefficients[[3]]; linear_slope <- summary(linear_lm)$coefficients[[2]]; linear_sse <- summary(linear_lm)$coefficients[[4]]     ushape_intercept <- summary(ushape_lm)$coefficients[[1]]; ushape_ise <- summary(ushape_lm)$coefficients[[3]]; ushape_slope <- summary(ushape_lm)$coefficients[[2]]; ushape_sse <- summary(ushape_lm)$coefficients[[4]]          z_intercept <- (ushape_intercept - linear_intercept) / sqrt(linear_ise^2 + ushape_ise^2)   #z = -0.45; p = .67, n.s.     z_slope <- (ushape_slope - linear_slope) / sqrt(linear_sse^2 + ushape_sse^2)   #z = -1.50; p = .93, sig.


title: Correlation using Logistic Regression and Pearson

text: I am so sorry, I am beginner in statistic analysis, I have project using R to analyze the correlation between dependent variables and independents variables. In this case I have two dependent variables (1. Extrovert, 2. Introvert). And the independent variables i have the data from (Call Log-> how long they call everyday, how many they call everyday, SMS log-> how length text in SMS body every day, how many they sent/received sms for each day). I am so confused how I can do it, please anyone can give me some good references about it. I also have some questions such as :   1. I use the different type of variables, independent variables (data type : numeric) but dependent variable (data type is categorical), so it is possible to apply logistic regression and Pearson?    2. Or any someone will give me some advice the better solution such as another methods for solving this problem.  The example of data from dput()               structure(list(sumcallin = c(462L, 998L, 335L, 179L, 34L, 0L,      0L, 0L, 0L, 0L), caountcallin = c(7L, 5L, 8L, 5L, 1L, 1L, 0L,      1L, 1L, 1L), sumcallout = c(1068L, 81L, 519L, 393L, 342L, 0L,      583L, 1902L, 358L, 1017L), countcallout = c(15L, 3L, 10L, 5L,      6L, 0L, 3L, 3L, 3L, 3L), sumreceived = c(322L, 75L, 20L, 35L,      8L, 35L, 135L, 103L, 471L, 173L), countreceived = c(15L, 4L,      2L, 3L, 1L, 2L, 7L, 3L, 18L, 5L), sumsent = c(171L, 31L, 25L,      23L, 8L, 55L, 87L, 9L, 400L, 258L), countsent = c(10L, 4L, 1L,      3L, 1L, 3L, 4L, 1L, 13L, 8L), personality = structure(c(2L, 2L,      2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("extro", "intro"), class = "factor")), .Names = c("sumcallin",      "caountcallin", "sumcallout", "countcallout", "sumreceived",      "countreceived", "sumsent", "countsent", "personality"), row.names = c(1L,      2L, 3L, 4L, 5L, 37L, 38L, 39L, 40L, 41L), class = "data.frame")      Thank you for your help.


title: Test if a slope falls within a back-transformed (log) prediction interval

text: I'm trying to test the hypothesis that the relationship (slope) between second molar tooth size and overall molar tooth size is 0.33 (in species of rodents), using generalized least squares regression in `R`. I'm doing this while "controlling" for autocorrelated residuals that occur because of the evolutionary relationships among species (using the `Gls` function in package `rms`). Both of my variables (`M2_area` and `total_area`) are quite right-skewed. When using arithmetic scale variables, the residuals from my model show distinct heterogeneity (`model1` below). I then refitted this model using inverse variance weights to try to improve this situation (`model2`), which it did to a limited degree. Since the variables are right-skewed and the residual variance increases with the size of the fitted values, I decided to try a natural log transformation for both predictor and response (`model3`). This seemed to me to produce the best diagnostic plots. But, i'm unsure how to test the above hypothesis using a log-log model. I know I can't back-transform the log-log coefficients. I have calculated predicted values and back-transformed these, together with their confidence intervals. A plot shows that the slope 0.33 does not lie within these confidence intervals. **My question: is there a way to formally test that the slope 0.33 does not lie inside the back-transformed prediction intervals?** (one caveat: I can't use any bootstrap procedures, as it is not possible to sample with replacement from an evolutionary tree, which is a necessary part of the GLS model) I include some data and R code below (I omit the correlation structure for the sake of brevity). # data               dat <- structure(list(M2_area = c(3.52, 3.66, 3.3, 2.65, 2.83, 3.39,      2.94, 2.74, 3.65, 3.02, 3.55, 3.46, 3.01, 2.97, 2.74, 3.69, 4.87,      1.46, 5.12, 4.61, 4.25, 5.28, 4.52, 3.41, 3.11, 2.81, 4.18, 2.54,      3.11, 2.58, 4.35, 6.59, 5.75, 4.21, 4.59, 3.82, 3.73, 3.64, 4.37,      3.79, 5.89, 6.31, 4.6, 4.18, 6.72, 5.63, 5.81, 5.19, 6.1, 5.58,      5.2, 4.82, 4.13, 6.1, 6.51, 12.05, 11.63, 12.94, 9.31, 12.16,      13.73, 11.35, 25.02, 11.24, 8.99, 12.24, 11.51, 16.07, 16.33,      4.92, 3.31, 2.98, 3.59, 3.43, 3.96, 6.07, 4.53, 4.21, 4.33, 4.96,      4.86, 2.55, 2.51, 4.55, 3.15, 0.35, 0.19, 0.69, 0.36, 0.36, 0.41,      0.9, 1.85, 1.75, 0.54, 0.56, 1.81, 1.08, 1.11, 1.19, 1.22), total_area = c(9.39,      10.17, 9.05, 7.23, 7.65, 9.38, 7.85, 7.49, 9.89, 8.16, 9.8, 9.33,      8.17, 7.88, 7.4, 9.98, 13.03, 3.95, 16.89, 14.95, 13.8, 17.52,      15.03, 11.28, 9.36, 8.15, 12.35, 8, 9.7, 8.54, 14.31, 20.65,      18.44, 13.81, 14.02, 12.43, 11.37, 11.16, 13.42, 12.32, 16.26,      18.84, 13.07, 12.78, 21.4, 17.11, 18.61, 16.42, 19.82, 17.69,      16.87, 14.48, 12.49, 20.15, 19.61, 36.74, 35.23, 40.31, 29.59,      37.64, 43.86, 37.17, 69.6, 33.99, 24.69, 33.71, 32.56, 45.63,      46.62, 13.44, 9.68, 8.13, 10.31, 10.43, 11, 16.97, 13.39, 12.96,      13.29, 14.27, 14.5, 7.27, 7.17, 12.89, 8.58, 0.86, 0.44, 1.55,      0.85, 0.88, 1.03, 2.54, 5.15, 5.16, 1.56, 1.61, 5.09, 3.01, 3.1,      3.42, 3.57)), .Names = c("M2_area", "total_area"), row.names = c(NA,      -101L), class = "data.frame")          # natural logs of variables     dat$ln_M2_area <- log(dat$M2_area)     dat$ln_total_area <- log(dat$total_area)      # models               library(rms)          # linear model (without corStruct for evolutionary tree)     model1 <- Gls(M2_area ~ total_area, data = dat)     model1     confint(model1)     qqnorm(model1)     plot(model1)          # linear model (with weighted least squares and without corStruct for tree)     model1fitted <- fitted(model1)     model2 <- Gls(M2_area ~ total_area, weights = varFixed(~ model1fitted),                   data = dat)     model2     confint(model2)     qqnorm(model2)     plot(model2)          # log-log model (without corStruct for tree)     model3 <- Gls(ln_M2_area ~ ln_total_area, data = dat)     model3     confint(model3)     qqnorm(model3)     plot(model3)          # create prediction interval for log-log model     newdat <- data.frame(ln_total_area = seq(min(dat$ln_total_area),                          max(dat$ln_total_area), by = 0.01))     pred_model3 <- data.frame(newdat, predict(model3, newdata = newdat,         type = "lp", se.fit = TRUE, conf.int = 0.95, conf.type = "mean"))     colnames(pred_model3)[2] <- "ln_M2_area"          # back transform logged prediction interval     smearing_estimate <- mean(exp(residuals(model3)))     pred_model3_back <- exp(pred_model3$ln_M2_area) * smearing_estimate     pred_model3_logci <- pred_model3$ln_M2_area +          (pred_model3$se.fit) %*% t(qnorm(c(0.5, 0.025, 0.975)))     pred_model3_ci <- exp(pred_model3_logci)     pred_model3_ci <- data.frame(exp(pred_model3$ln_total_area), pred_model3_ci)     colnames(pred_model3_ci) <- c("total_area", "M2_area", "lcl", "ucl")     head(pred_model3_ci)      # plot               library(ggplot2)          ggplot(dat, aes(x = total_area, y = M2_area)) +          geom_point(shape = 1, size = 2.2, alpha = 0.6) +         geom_line(data = pred_model3_ci, color = "red") +          geom_ribbon(data = pred_model3_ci, aes(ymin = lcl, ymax = ucl),                      alpha = 0.2) +         geom_abline(intercept = 0, slope = 0.333, linetype = "dashed",                      color = "grey30") +         coord_fixed(ratio = 1) +         theme_bw()


title: Extracting metadata (author, creation date ...) from word documents with R

text: I would like to extract the metadata (author, creation date ...) from a batch of MS word documents with R. Can anybody tell me if that is possibile, and if so how I can do this? Thanks, Fred.


title: Binomial data and PCA and cluster analysis

text: I have obtained responses on around 48 items measuring employer attractiveness. The goal of my analysis is to cluster respondents according to these employer attractiveness dimensions. I plan to reduce the dimensionality with PCA and also confirm the structure with CFA before getting to cluster analysis. 48 items were measured in forced 4-point Likert scale, but the data I obtained for analysis was transformed by another research to binomial level - important or not important (0 or 1 values). I was wondering what effect that can have on PCA and cluster analysis that these data are binomial. Can these analysis techniques be executed on binomial data?


title: Question about basic moderation analysis

text: Assume a 2x2 between-subjects experimental design (FACTOR1 and FACTOR2, which are both categorical), and there's also 2 continuous dependent variables, DV1 and DV2. Assume that on ANOVA on DV1 reveals that there's no sig. effect except a main effect of FACTOR1, whereas the ANOVA on DV2 reveals no sig. effect. What I'm wondering is, would it be appropriate to do a moderation analysis IN ADDITION to these ANOVAs? The moderation analysis is to check whether the relation between DV1 and DV2 is moderated by FACTOR1.


title: How to calculate survey sample with more than one outcome variables

text: I'm working on a proposal and that needs survey design. The basic question is like this: The study seeks to examine factors (e.g., income) associated with the use of health checkup services of workers at an enterprise (they have to pay out of pocket for these services). Among those who check up, some will have health issues and others will not. The study also wants to examine factors influencing use of further, advanced services and treatments in the group encountering health issues. So there are two outcomes, and one population for analysis (checkup workers) is a subset of a bigger population (workers of the enterprise). I'm quite new to sampling methods, so I wonder if you can give some useful readings for sampling strategy in this situation. I mean, many sampling textbooks are out there, but I want to know what topics and key words I need to focus on and read.


title: Estimation/identification simple problem

text: I have this model in my mind: $ y_t = c + \phi x_{t-1}^k + \eta_t$ where $\eta_t$ is N.I.D. In your opinion is possible to estimate $k$?


title: How to calculate percent variation

text: I've came across a method in research communication, where different percentages of variation were calculated for a collection of ratio values for some proteins. As seen in the figure, the frequency of (protein) ratios falling within a certain percentage of variation is plotted against the percent variation itself. ![http://i.stack.imgur.com/RgXkC.jpg](http://i.stack.imgur.com/5u3DR.jpg) Can anyone tell me how this percent variation is calculated?


title: Calculating precision and recall in R

text: Suppose I'm building a logistic regression classifier that predicts whether someone is married or single. (1 = married, 0 = single) I want to choose a point on the precision-recall curve that gives me at least 75% precision, so I want to choose thresholds $t_1$ and $t_2$, so that:   * If the output of my classifier is greater than $t_1$, I output "married".   * If the output is below $t_2$, I output "single".   * If the output is in between, I output "I don't know". A couple questions:   1. I think under the standard definition of precision, precision will be measuring the precision of the married class alone (i.e., precision = # times I correctly predict married / total # times I predict married). However, what I really want to do is measure the overall precision (i.e., the total # times I correctly predict married or single / total # times I predict married or single). Is this an okay thing to do? If not, what should I be doing?   2. Is there a way to calculate this "overall" precision/recall curve in R (e.g., using the ROCR package or some other library)? I'm currently using the ROCR package, but it seems to only give me the single-class-at-a-time precision/recall.


title: Is a multigroup confirmatory factor analysis appropriate for comparing the measurement model across two groups?

text: **Background:** Before commencing treatment 180 participants completed a baseline 17-item Hamilton depression rating scale (HDRS), which is a likert scale. Because of the treatment, half developed depression and the other half did not. So i have two groups; depressed and non-depressed. **Analysis:** I want to know whether the factor structure of baseline HDRS differs between the two groups. Thus, I performed multiple group confirmatory factor analysis (CFA). I used a theoretically validated and relevant HDRS factor model as the basis of the MGCFA to compare goodness of fit measures and determine measurement invariances between the groups. However, my initial CFA indicated a poor model fit for both my groups. I am mindful that my sample sizes may be too small for this. ### Questions   * Is multigroup CFA appropriate for comparing the measurement model of my two groups?   * What should I do given the poor model fit of the CFA in both groups?   * Is a total sample size of 180 too small for this multiple group CFA?


title: Alternatives to classification trees, with better predictive (e.g: CV) performance?

text: I am looking for an alternative to Classification Trees which might yield better predictive power. The data I am dealing with has factors for both the explanatory and the explained variables. I remember coming across random forests and neural networks in this context, although never tried them before, are there another good candidate for such a modeling task (in R, obviously)?


title: Comparing Unemployment Time Series

text: I am looking to see if recessions affect the unemployment rate in capital cities differently than the state as a whole. I have the data for the unemployment rate each month from January 2007 to July 2008 of all 50 states, their capitals, and the capital metropolitan area. I primarily am looking to see if one lags behind another (ex. capital unemployment rate on average drops slower than the states), or if one has more drastic changes than another (ex. state's unemployment rate on average drops more than capitals). As you can probably tell by my vocabulary, I am a complete amateur at statistics, but any help in pointing me in the right direction would be appreciated. Also, could such analysis be done using Excel? Thanks in advance. What my data looks like if I confused anyone:                                   Jan-07  Feb-07  Mar-07  Apr-07     Arizona              3.8    3.7     3.7      3.5     Arizona Capital      4.3    3.8     3.6      3.4     Arizona Capital M    3.7    3.3     3.1      2.9     Arkansas             5.2    5.2     5.2      5.3     Arkansas Capital     5.3    5.4     4.9      4.8     Arkansas Capital M   4.9    5.0     4.5      4.4


title: Being totally stuck at my master thesis: cant figure out this 3 equation SEM

text: I got smth which looks like this $y_1 = \alpha_1\cdot y_2 + \alpha_2\cdot y_3 + X\cdot\alpha_3 + u_1$ $y_2 = \beta_1\cdot y_1 + \beta_2\cdot y_3 + X\cdot\beta_3 + u_2$ $y_3 = \gamma_1\cdot y_1 + \gamma_2\cdot y_2 + X\cdot\gamma_3 + u_3$ After I put some thought into my model I could not really restrict some of the exogenous variables (which are contained in $X$) to zero nor could I kick out some of the dependencies between those endogenous variables... The topic which my model is about will regard to the success of a student during his studies. Hence the dependent variables will be grade, semester and how many times a student revisited a exam (because of fail or to get a higher grade). First I tried IV, 2SLS, 3SLS but since my data just do not fit the requirements of this methods this is a dead end. The reason is that a) there is no a good instrument available and b) even if I try to find the reasons to not include one of the exogenous variables in one or more equations the first stage fit is so bad that the second stages suffers a great deal from it. The next thing I tried was to introduce a approach which allows me to not exclude any variables by including a rank-order variable (introduced by vella). This RO-Variable is supposed to account for the endogeneity by transforming the reduced form residuals such that the new error term in the structural equation is now uncorrelated with the endogenous regressors. The reason this did not work out was that the costs of this approach lies in the high multicolinearity which renders my results useless. After I did a simulation I figured out that my results pointed towards some heavy problems which cant be overcome. I think this is a very neat solution which was the reason I took it but it just cannot work out. So the bottom line is: I'm stuck, I don't now what I could do next and I hope for some input of you guys what I could possibly do to make this work. I just need some hints in the right direction because I feel a bit lost..


title: Adding not chosen alternatives as data to logistic regression model

text: I am interested in predicting shopping behaviour in a shopping center. I have a database with the chosen alternatives (shop) and variables describing that alternative (like type and size) and the individual which made the choice (sex, age). What I would like to get is the influcence of these variables on the behaviour. My question: Would it be possible to append for each chosen shop the not chosen alternatives and calculate the binomial logistic model? As an example:               | person.id | shop.id | chosen | type    | size | sex | age |     |         1 | A       | yes    | fashion |    3 | m   |  21 |     |         1 | B       | no     | grocery |    5 | m   |  21 |     |         1 | C       | no     | fashion |    1 | m   |  21 |     |         2 | B       | yes    | grocery |    5 | f   |  45 |     |         2 | A       | no     | fashion |    3 | f   |  45 |     |         2 | C       | no     | fasion  |    1 | f   |  45 |      I was first thinking about the multinomial logistic regression but I am actually only interested in the variable chosen which is dichotomous (yes/no).


title: nth moment, for 0 < n < 1 or n <0, do they exist?

text: I am interested in the moments, we have for instance the mean, $\mathrm{E}(X)$ and $\mathrm{E}(X^2)$. What about values like $\mathrm{E}(X^{1.5})$ or $\mathrm{E}(X^{-1})$? Have they been investigated?


title: directed bayesian network and factor graphs

text: I have a directed bayesian given by the figure below. In the figure the circles are random variables and the shaded ones are observed. The rectangular nodes are constants representing the hyper parameters over the prior distributions. ![Probabilistic graphical model](http://i.stack.imgur.com/9Spyz.png) I have now tried to prepare the corresponding factor graph for the problem, which I hope represents the same joint distribution as the product of the factors (denoted by the rectangular nodes). The random variables are denoted similarly using circles. This is shown in the image below. Can someone comment on whether I have done it correctly? ![Corresponding factor graph](http://i.stack.imgur.com/yQxDb.png) So, as far as I think, I have the following: $$ f_1 = P(\lambda) \\\ f_2 = P(w|\lambda) \\\ f_3 = P(y|w, \phi, x) \\\ f_4 = P(\phi) $$ So, the product represents $P(\lambda) \times P(w|\lambda) \times P(y|w, \phi, x) \times P(\phi)$, which is the joint distribution specified by the original factor graph i.e. $$ f(\lambda, \phi, w) = f_1(\lambda) f_2(w, \lambda) f_3(w, \phi) f_4(\phi) $$ Also, now if I want to compute the marginal at $w$, can I write it as: $$ P(w) = \sum_{\lambda}f_2(w, \lambda)f_1(\lambda) \sum_{phi}f_3(w, \phi) f_4(\phi) $$ The idea being I compute this joint distribution and then I can marginalise by summing and product of factors somehow.


title: how to interpret KS test or Shapiro wilk test for ordinal criterion variable?

text: ![enter image description here](http://i.stack.imgur.com/0PBg0.png) Here Rank is my dependent variable from level (Strongly disagree to strongly agree) and Pre_Sales_support is one of the independent variable. Due to smaller sample size of 15, I am taking shapiro-wilk test into account to check normality but my null hypothesis getting rejected for one level(p value for Strongly agreed .007) and accepted for another(p value for agreed .101) within same variable. How to get around this.


title: Calculating standard error and attaching an error bar on ggplot2 bar chart

text: Given a minimal dataset where am looking for the occurrence of a certain motif within a dataset of 500 observations. with_motif represents obervations with the specified motif and without_motif are observations without the motif.               with_motif <- 100     without_motif <- 400     dt <- data.frame(with_motif,without_motif)      The following code will plot a bar-chart using ggplot2 library,               bar_plot <- ggplot(melt(dt),aes(variable,value)) + geom_bar() + scale_x_discrete(name="with or without") + theme_bw() + opts( panel.grid.major = theme_blank(),title = "", plot.title=theme_text(size=14))          bar_plot      I would like to compute a standard error at 95% CI and attach a barchart to the plot. ggplot offers `geom_errorbar()` but I would be glad to know different ways for deriving the standard errors(deviation) so as to calculate the errorbar limits(CI).


title: Why are degrees of freedom so high in repeated measures mixed models?

text: I have been running a repeated measures mixed model in SAS and the denominator degrees of freedom for the F statistic seems to be really high. In this case I have roughly 200 individuals and three measurements for each individual (the repeated statement) and degrees of freedom are over 500. So it seems to be roughly reflecting the number of observations. I checked a few articles and found the same pattern elsewhere, but as a stats beginner I'm not sure if it is ok; or is it pseudoreplication? I also noted that changing the covariate structure has a large effect on degrees of freedom.   * Why are degrees of freedom so high in such designs, or is it likely that there is something seriously wrong with my model?   * What might be some suggested reading on repeated measures mixed models to help me understand this?


title: Real based on Power Function

text: Problem:   What is a real life example of a power function? I tried thinking about it but I failed. Does somebody know?


title: Nonparameteric multivariate density approximation -- where do I start?

text: I am currently working on a research project that requires a reliable method for non-parametric kernel density estimation. Some specifics about my problem:   * I have $N$ sample points $X_1,X_2...X_N$, each of which corresponds to an outcome of an $M$ dimensional random variable $W = (w_1..w_M)$. I would like to assign a density to each of these points in an accurate and consistent way.    * I am looking for a method that does not assume that the random variables be independent. As in my situation, I know that $f(W) \neq \prod_{i=1}^M {f_m(w_i)}$.   * In my case, $N \leq 25000$ and $M \leq 5$ but I can generate as much data as needed.  Open to any and all suggestions!


title: Law of iterated expectations - an small exercise

text: From edX MIT probability course - Widgets and Crates: * * * Let $X_i$ be the number of widgets in a particular box $i$.   Let $N$ be the number of boxes in a crate. Assume $X$ and $N$ are independent, both with expected value equal to 10 and variance equal to 16. Evaluate the expected value of $T$, where $T$ is the total widgets in a crate. * * * Is it wrong to take the following aproach? (1) Law of iterated expectations: $$E[T] = E(E[T|N])$$ (2) Since in a crate there will be $X.N$ widgets, with both $X$ and $N$ random, $$E[E(T|N)]= E[N.X]$$ (3) Since $X$ and $N$ are independent, $$E[X.N] = E[X].E[N]=10^2=100$$       Thanks in advance,   SB


title: Suspicious Friedman test results

text: My data involves testing 216 experimental units (computational models) for which one factor (the value of a specific shared parameter of the models) is administered in 4 treatments (i.e. the parameter possible values). The experimental design is repeated measures i.e. each model is tested with all the 4 values of the parameters. My observations corresponds to the "accuracy" of the model given a certain parameter value (or actually the AUC of the ROC analysis curve). There are no missing values, however there are many ties (i.e. many outcomes coincide). My hypothesis is that the value of this parameter does not influence significantly the AUC of the model; or in other words, that the variation in the model's AUC given the parameter value is attributable to chance. In a sense, I will be expecting in a non-significant result of my hypothesis testing. ![Boxplots](http://i.stack.imgur.com/pYaip.png) Now, following visual exploration of my data, and representation of the data distribution factorized by the treatment in the form of notched boxplots, intuitively any test would answer negatively (considerable overlap of the factorized distributions, especially on the notched region proportional to the standard error). However, when running Friedman in R I got the following results: Friedman chi-squared = 50.3378, df = 3, p-value = 6.77e-11 Suspicious of this result, I tested in a tool I found online (http://vassarstats.net/index.html), and the results were: Friedman chi-squared = 13.94, df = 3, p-value = 0.003 I found related information about this disagreement in here: http://r.789695.n4.nabble.com/Unexpected-behavior-in-friedman-test-and-ks- test-td902324.html Still, based on the boxplots, I'm unconvinced that any test shall result in a significant finding. Perhaps the many ties are be fooling the implementations?. Also, I may be overpowering with too many observations, thus will obtain a significant result regardless of the data. The test assumptions seem to hold (unless I'm totally wrong). A bug in the R code seems unlikely after so many checkings (but can never be discarded). My question is; in the light of the boxplots; why does Friedman gives a significant result? or am I doing something wrong? or am I being silly in applying Friedman when I should be applying other test? or should I trust this finding (blindly)? I would appreciate plain and direct answers as I'm not an statistician.


title: Let $f(x)$ be the density of $X$. What is $\lim_{n\rightarrow \infty}\mathrm{E}_X[nf(X+n)]$?

text: The expectation can be written as $\mathrm{E}_X[nf(X+n)] = \int_{-\infty}^\infty nf(x+n)f(x)\,\mathrm{d}x$. The expectation relies on the speed of tail $f(x+n)$ goes to $0$. I posted a related question in this link.


title: Regression line fit for linearly increasing data with manual reset

text: I've a linearly increasing time series dataset of a sensor with value ranges between 50 to 150 on which I've implemented a simple linear regression algorithm to fit a regression line, and I'm predicting the date when the series would reach 120, 130 and 140. All works fine when the series move upwards, but there are cases when the sensor reaches around 110 or 115 in which a person would manually reset the sensor, and the values start over again at say 50 or 60... This is where I start facing issues with the regression line as it starts moving downwards, and it starts predicting old date. I think I should be considering only the subset of data from where it was previously reset. However, I'm trying to understand if there are any algorithms available that consider this case. I'm new to data science, and would appreciate any pointers to move further. Thanks.


title: R: How to to simulate ARIMA using starting values?

text: I have built an ARIMA(p,d,q) model, m using say,               m <- Arima(ts.data, c(p,d,q))      Given some starting values, I want to simulate future values based on the model m. I can manually code it but there must be a way to simulate it. In `arima.sim` there doesn't seem to be a way to specify the starting values.


title: What is the right way to analyze a nested design in R?

text: I know that there are already a host of questions about nested designs but many of them haven't been answered or come from biological domains which I sometimes find hard to transfer to my domain. I am currently trying to analyse the data from a psychological study. Participants read four statements and indicated how appropriate they found each on a 5-point Likert-scale. But I'll just call the outcome variable "out" and so the content will not distract from the statistical problem at hand. There are two IVs, one that was varied on three levels between subjects (let's call it "between" and the levels 1, 2 and 3) and one on two levels within subjects ("within" with its levels a and b). The within subjects factor is realised in such a manner that two of the four statements correspond to one of the levels and the other two to the other. Correspondingly, for every participant, I now have four "out" data points, two for the a statements and two for the b statements, like this (in long format):               subj    between within  item    out     123     1       a       a1      3     123     1       a       a2      4     123     1       b       b1      1     123     1       b       b2      2     124     2       a       a1      5     124     2       a       a2      4     124     2       b       b1      2     124     2       b       b2      3     125     3       a       a1      1     125     3       a       a2      1     125     3       b       b1      2     125     3       b       b2      3      and so on What is the right way to analyse the effect and interaction of between and within? I assume that the two different a and the two different b statements do not differ systematically but are just two instances of the same thing. Or do I need an additional factor that tells R which of the four statements it is in order to allow for that error? Anyway, I have tried these options:               m1 <- aov(out~between*within+Error(subj/within))     summary(m1)          m2 <- lm(out~between*within+within/subj)     anova(m2)          m3 <- lmer(out~between*within+(between|subj))     Anova(m3)      They produce different results but I'm not entirely sure about which one is right or what the differences are. Could anyone enlighten me on this? I assume this is about fixed and random effects which always gets my head in a twist. I have read other posts about nested designs and I think this one comes nearest. But unfortunately the question ID being nested groups wasn't answered. Any help would be very much appreciated! **Edit:**   I greatly appreciate the answers so far! The comments prompted me to read some more tutorials on mixed effects models as well as answers to similar questions. This enables me to clarify my question: I know that I will need to specify random intercepts at least for subjects and probably also the four different items representing the two levels of the within factors. However, I am unsure about the specific fixed/random effects structure I would have to model, and specifically what the maximal model in my case would look like. Right now, I have tried some code which roughly looks like this:               m1 <- lmer(out ~ between * within + (between|subj) + (within|subj) + (between|within/item))      This is probably wrong on several levels so any feedback would be greatly appreciated! Just for clarity: The item factor has no intrinsic meaning, but both levels of the within factor is realised by two items each (which every participant sees) and I am unsure whether this is relevant to the model. Apart from the model specification, I am still getting confused about the concepts nested and crossed and which one applies to my setup. Lastly, most of my models fail to converge but that is probably a matter for another question (and has been discussed here in length.)


title: Convergence of Estimation

text: The issue is, we can have a bunch of regression models to do prediction for the test set. For a specific testing data, some perform good and some do bad. How can I use existing method or design a convergence precess so that **the prediction could converge to the "best" response after using a part of models**? For example, we have many sub decision trees when applying random forest. In random forest, all sub trees are used to do prediction, then voted for the best response. If we can have convergence process, maybe we just need some of them and converge to the best response in the end. Help~